Cod diet: analysis of Levins’s diversity index and Schoener’s overlap index

Intro

Calculate Schoener’s overlap and Levin’s diversity index and fit brms models

# Load libraries, install if needed
library(tidyverse); theme_set(theme_light(base_size = 12))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(EnvStats)
library(qwraps2)
library(bayesplot)
library(tidybayes)
library(brms)

# To load entire cache in interactive r session, do: qwraps2::lazyload_cache_dir(path = "R/analysis/cod_flounder_diets_spatial_cache/html")

For maps

# Specify map ranges
ymin = 55; ymax = 58; xmin = 14; xmax = 20

map_data <- rnaturalearth::ne_countries(
  scale = "medium",
  returnclass = "sf", continent = "europe")

# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
  st_crop(map_data,
          c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))

# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)

ggplot(swe_coast_proj) + geom_sf() 

# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
      axis.text.x = element_text(angle = 90),
      axis.text = element_text(size = 8),
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 8),
      legend.position = "bottom",
      legend.key.height = unit(0.2, "cm"),
      legend.margin = margin(0, 0, 0, 0),
      legend.box.margin = margin(-5, -5, -5, -5),
      strip.text = element_text(size = 8, colour = 'black', margin = margin()),
      strip.background = element_rect(fill = "grey90")
      )
}

# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
        axis.text.x = element_text(angle = 90),
        axis.text = element_text(size = 6),
        strip.text = element_text(size = 8, colour = 'black', margin = margin()),
        strip.background = element_rect(fill = "grey90")
      )
}

# Make default base map plot
plot_map_raster <- 
ggplot(swe_coast_proj) + 
  geom_sf(size = 0.3) +
  labs(x = "Longitude", y = "Latitude") +
  theme_facet_map(base_size = 14)

# Function to convert from lat lon to UTM
LongLatToUTM <- function(x, y, zone){
  xy <- data.frame(ID = 1:length(x), X = x, Y = y)
  coordinates(xy) <- c("X", "Y")
  proj4string(xy) <- CRS("+proj=longlat +datum=WGS84")  ## for example
  res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
  return(as.data.frame(res))
}

theme_clean <- function() {
  theme_minimal(base_family = "Barlow Semi Condensed") +
    theme(panel.grid.minor = element_blank(),
          plot.title = element_text(family = "BarlowSemiCondensed-Bold"),
          axis.title = element_text(family = "BarlowSemiCondensed-Medium"),
          strip.text = element_text(family = "BarlowSemiCondensed-Bold",
                                    size = rel(1), hjust = 0),
          strip.background = element_rect(fill = "grey80", color = NA))
}

Read and plot data

# This is just to add the density covariates to the schoener overlap models
cod <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/cod_diet_analysis.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#> 
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#>   .default = col_double(),
#>   unique_pred_id = col_character(),
#>   time_period = col_character(),
#>   quarter_fact = col_character(),
#>   pred_weight_source = col_character(),
#>   predator_spec = col_character(),
#>   ices_rect = col_character(),
#>   cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.

cod <- cod %>%
  mutate(year = as.integer(year),
         quarter = as.factor(quarter),
         depth2_sc = depth - mean(depth),
         saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
         tot_prey_biom_per_mass = tot_prey_biom/pred_weight_g,
         depth_sc = (depth - mean(depth)) / sd(depth)) %>% 
  filter(year > 2014) %>%
  filter(!quarter == 2) %>% 
  drop_na(predfle_density_sc, predcod_density_sc) %>% 
  droplevels()
#> mutate: converted 'year' from double to integer (0 new NA)
#>         converted 'quarter' from double to factor (0 new NA)
#>         new variable 'depth2_sc' (double) with 102 unique values and 0% NA
#>         new variable 'saduria_entomon_per_mass' (double) with 1,666 unique values and 0% NA
#>         new variable 'tot_prey_biom_per_mass' (double) with 7,108 unique values and 0% NA
#>         new variable 'depth_sc' (double) with 102 unique values and 0% NA
#> filter: removed 4,936 rows (58%), 3,593 rows remaining
#> filter: no rows removed
#> drop_na: no rows removed

# Now read stomach data (1 row 1 stomach)
cod_prey <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/cod_diet_full_analysis.csv") %>% dplyr::select(-X1) %>% mutate(species = "COD") %>% filter(year > 2014) %>% filter(!quarter == 2) %>% droplevels()
#> Warning: Missing column names filled in: 'X1' [1]
#> 
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#>   .default = col_double(),
#>   unique_pred_id = col_character(),
#>   predator_spec = col_character(),
#>   ices_rect = col_character(),
#>   cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> filter: removed 4,936 rows (58%), 3,593 rows remaining
#> filter: no rows removed

fle_prey <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/fle_diet_full_analysis.csv") %>% dplyr::select(-X1) %>% mutate(species = "FLE") %>% filter(!quarter == 2) %>% droplevels()
#> Warning: Missing column names filled in: 'X1' [1]
#> 
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#>   .default = col_double(),
#>   unique_pred_id = col_character(),
#>   predator_spec = col_character(),
#>   ices_rect = col_character(),
#>   cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> filter: no rows removed

# Plot data in space
plot_map_raster +
  geom_point(data = fle_prey, aes(x = X * 1000, y = Y * 1000, color = "FLE"), alpha = 0.5) +
  geom_point(data = cod_prey, aes(x = X * 1000, y = Y * 1000, color = "COD"), alpha = 0.5) +
  theme_plot() +
  facet_wrap(~sub_div)

cod_prey <- cod_prey %>% filter(lat < 58)
#> filter: no rows removed
fle_prey <- fle_prey %>% filter(lat < 58)
#> filter: no rows removed

# I will also make a new area category, pooling 24 and 25, and 27-8, making it a southwest (coastal), north (coastal) and offshore (26). Note this is only for the cod_prey and fle_prey data, which we use for indices. The other stomach analysis is spatial, but here we need enough samples to calculate indices properly
cod_prey <- cod_prey %>% mutate(area = ifelse(sub_div %in% c(24, 25), "24-25", NA),
                                area = ifelse(sub_div %in% c(27, 28), "27, 28", area),
                                area = ifelse(sub_div == 26, "26", area))
#> mutate: new variable 'area' (character) with 3 unique values and 0% NA

fle_prey <- fle_prey %>% mutate(area = ifelse(sub_div %in% c(24, 25), "24-25", NA),
                                area = ifelse(sub_div %in% c(27, 28), "27, 28", area),
                                area = ifelse(sub_div == 26, "26", area))
#> mutate: new variable 'area' (character) with 3 unique values and 0% NA

Calculate Levin’s index for diet (by species, all sizes pooled and by area to ensure large enough sample size)

# Reformat data
# colnames(fle_prey)
fle_prey_long <- fle_prey %>%
  pivot_longer(15:30) %>% 
  group_by(name, year, quarter, area) %>%
  summarise(fle_stomach_content = sum(value)) %>% 
  arrange(name, year, area, quarter) %>%
  ungroup() %>% 
  mutate(id = paste(year, quarter, area, sep = "_")) %>% 
  group_by(id) %>%
  mutate(fle_stomach_content_tot = sum(fle_stomach_content),
         fle_stomach_content_prop = fle_stomach_content / fle_stomach_content_tot) %>% 
  ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 2180x32, now 34880x18]
#> group_by: 4 grouping variables (name, year, quarter, area)
#> summarise: now 352 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 22 unique values and 0% NA
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'fle_stomach_content_tot' (double) with 22 unique values and 0% NA
#>                   new variable 'fle_stomach_content_prop' (double) with 199 unique values and 0% NA
#> ungroup: no grouping variables

# This should amount to the number of unique prey
number_of_prey <- fle_prey_long %>% group_by(id) %>% mutate(n = n()) %>% ungroup() %>% distinct(n)
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'n' (integer) with one unique value and 0% NA
#> ungroup: no grouping variables
#> distinct: removed 351 rows (>99%), one row remaining
number_of_prey
#> # A tibble: 1 × 1
#>       n
#>   <int>
#> 1    16

test_id <- head(fle_prey_long$id, 1)
fle_prey_long %>% filter(id == test_id) %>% as.data.frame()
#> filter: removed 336 rows (95%), 16 rows remaining
#>                      name year quarter  area fle_stomach_content           id
#> 1           amphipoda_tot 2015       4 24-25                7.64 2015_4_24-25
#> 2            bivalvia_tot 2015       4 24-25              291.27 2015_4_24-25
#> 3     clupea_harengus_tot 2015       4 24-25                0.00 2015_4_24-25
#> 4           clupeidae_tot 2015       4 24-25                0.00 2015_4_24-25
#> 5          gadiformes_tot 2015       4 24-25                0.00 2015_4_24-25
#> 6        gadus_morhua_tot 2015       4 24-25                0.00 2015_4_24-25
#> 7            gobiidae_tot 2015       4 24-25                0.00 2015_4_24-25
#> 8             mysidae_tot 2015       4 24-25                0.00 2015_4_24-25
#> 9             non_bio_tot 2015       4 24-25                6.59 2015_4_24-25
#> 10    other_crustacea_tot 2015       4 24-25                4.05 2015_4_24-25
#> 11       other_pisces_tot 2015       4 24-25                0.00 2015_4_24-25
#> 12              other_tot 2015       4 24-25                3.44 2015_4_24-25
#> 13 platichthys_flesus_tot 2015       4 24-25                0.00 2015_4_24-25
#> 14         polychaeta_tot 2015       4 24-25                2.82 2015_4_24-25
#> 15    saduria_entomon_tot 2015       4 24-25               45.85 2015_4_24-25
#> 16  sprattus_sprattus_tot 2015       4 24-25                0.00 2015_4_24-25
#>    fle_stomach_content_tot fle_stomach_content_prop
#> 1                   361.66              0.021124813
#> 2                   361.66              0.805369684
#> 3                   361.66              0.000000000
#> 4                   361.66              0.000000000
#> 5                   361.66              0.000000000
#> 6                   361.66              0.000000000
#> 7                   361.66              0.000000000
#> 8                   361.66              0.000000000
#> 9                   361.66              0.018221534
#> 10                  361.66              0.011198363
#> 11                  361.66              0.000000000
#> 12                  361.66              0.009511696
#> 13                  361.66              0.000000000
#> 14                  361.66              0.007797379
#> 15                  361.66              0.126776530
#> 16                  361.66              0.000000000

# Now cod
# colnames(cod_prey)
cod_prey_long <- cod_prey %>%
  pivot_longer(15:30) %>% 
  group_by(name, year, quarter, area) %>%
  summarise(cod_stomach_content = sum(value)) %>% 
  arrange(name, year, area, quarter) %>%
  ungroup() %>% 
  mutate(id = paste(year, quarter, area, sep = "_")) %>% 
  group_by(id) %>%
  mutate(cod_stomach_content_tot = sum(cod_stomach_content),
         cod_stomach_content_prop = cod_stomach_content / cod_stomach_content_tot) %>% 
  ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x32, now 57488x18]
#> group_by: 4 grouping variables (name, year, quarter, area)
#> summarise: now 416 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 26 unique values and 0% NA
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'cod_stomach_content_tot' (double) with 26 unique values and 0% NA
#>                   new variable 'cod_stomach_content_prop' (double) with 299 unique values and 0% NA
#> ungroup: no grouping variables

unique(is.na(fle_prey_long))
#>       name  year quarter  area fle_stomach_content    id
#> [1,] FALSE FALSE   FALSE FALSE               FALSE FALSE
#>      fle_stomach_content_tot fle_stomach_content_prop
#> [1,]                   FALSE                    FALSE
unique(is.na(cod_prey_long))
#>       name  year quarter  area cod_stomach_content    id
#> [1,] FALSE FALSE   FALSE FALSE               FALSE FALSE
#>      cod_stomach_content_tot cod_stomach_content_prop
#> [1,]                   FALSE                    FALSE

fle_prey_long %>% mutate(fle_stomach_content_prop = replace_na(fle_stomach_content_prop, -9)) %>% filter(fle_stomach_content_prop == -9)
#> mutate: no changes
#> filter: removed all rows (100%)
#> # A tibble: 0 × 8
#> # … with 8 variables: name <chr>, year <dbl>, quarter <dbl>, area <chr>,
#> #   fle_stomach_content <dbl>, id <chr>, fle_stomach_content_tot <dbl>,
#> #   fle_stomach_content_prop <dbl>
fle_prey_long <- fle_prey_long %>% mutate(fle_stomach_content_prop = replace_na(fle_stomach_content_prop, 0)) 
#> mutate: no changes

# Calculate Levin niche index
levin <- left_join(cod_prey_long, fle_prey_long) %>% 
  drop_na(name) %>% 
  drop_na(fle_stomach_content_prop) %>% 
  drop_na(cod_stomach_content_prop) %>% 
  group_by(year, quarter, area) %>% 
  summarise(levin_cod = (1/(number_of_prey$n-1)) * (((1/sum(cod_stomach_content_prop^2))) - 1),
            levin_fle = (1/(number_of_prey$n-1)) * (((1/sum(fle_stomach_content_prop^2))) - 1)) %>% 
  ungroup()
#> Joining, by = c("name", "year", "quarter", "area", "id")
#> left_join: added 3 columns (fle_stomach_content, fle_stomach_content_tot, fle_stomach_content_prop)
#>            > rows only in x    64
#>            > rows only in y  (  0)
#>            > matched rows     352
#>            >                 =====
#>            > rows total       416
#> drop_na: no rows removed
#> drop_na: removed 64 rows (15%), 352 rows remaining
#> drop_na: no rows removed
#> group_by: 3 grouping variables (year, quarter, area)
#> summarise: now 22 rows and 5 columns, 2 group variables remaining (year, quarter)
#> ungroup: no grouping variables

levin
#> # A tibble: 22 × 5
#>     year quarter area   levin_cod levin_fle
#>    <dbl>   <dbl> <chr>      <dbl>     <dbl>
#>  1  2015       4 24-25     0.216     0.0335
#>  2  2015       4 26        0.111     0.152 
#>  3  2015       4 27, 28    0.232     0.0623
#>  4  2016       1 24-25     0.154     0.174 
#>  5  2016       1 26        0.165     0.0352
#>  6  2016       1 27, 28    0.0736    0.112 
#>  7  2016       4 24-25     0.112     0.0231
#>  8  2016       4 26        0.137     0.0292
#>  9  2016       4 27, 28    0.0780    0.122 
#> 10  2017       1 24-25     0.146     0.0602
#> # … with 12 more rows

# Quickly check data to determine which distribution to use
levin %>% 
  ungroup() %>% 
  count(levin_cod == 0) %>% 
  mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now one row and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with one unique value and 0% NA
#> # A tibble: 1 × 3
#>   `levin_cod == 0`     n  prop
#>   <lgl>            <int> <dbl>
#> 1 FALSE               22     1

levin %>% 
  ungroup() %>% 
  count(levin_fle == 0) %>% 
  mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now one row and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with one unique value and 0% NA
#> # A tibble: 1 × 3
#>   `levin_fle == 0`     n  prop
#>   <lgl>            <int> <dbl>
#> 1 FALSE               22     1

# Fit beta models, so few zeroes and no 1's
levin %>% arrange(desc(levin_cod)) %>% dplyr::select(levin_cod)
#> # A tibble: 22 × 1
#>    levin_cod
#>        <dbl>
#>  1     0.246
#>  2     0.232
#>  3     0.216
#>  4     0.188
#>  5     0.165
#>  6     0.159
#>  7     0.154
#>  8     0.149
#>  9     0.146
#> 10     0.137
#> # … with 12 more rows
levin %>% arrange(desc(levin_fle)) %>% dplyr::select(levin_fle)
#> # A tibble: 22 × 1
#>    levin_fle
#>        <dbl>
#>  1    0.184 
#>  2    0.174 
#>  3    0.166 
#>  4    0.153 
#>  5    0.152 
#>  6    0.129 
#>  7    0.122 
#>  8    0.114 
#>  9    0.112 
#> 10    0.0993
#> # … with 12 more rows

# Final polish of data before feeding into models (species, not size-based indicies)
levin <- levin %>%
  mutate(levin_cod2 = ifelse(levin_cod == 0, 0.0001, levin_cod),
         levin_fle2 = ifelse(levin_fle == 0, 0.0001, levin_fle),
         year_f = as.factor(year),
         quarter_f = as.factor(quarter),
         area_f = as.factor(area))
#> mutate: new variable 'levin_cod2' (double) with 22 unique values and 0% NA
#>         new variable 'levin_fle2' (double) with 22 unique values and 0% NA
#>         new variable 'year_f' (factor) with 6 unique values and 0% NA
#>         new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#>         new variable 'area_f' (factor) with 3 unique values and 0% NA

Calculate Schoener’s index (by species)

# Reformat data to calculate Schoeners overlap index
# colnames(fle_prey)
fle_prey_long <- fle_prey %>%
  pivot_longer(15:30) %>% 
  group_by(name, year, quarter, ices_rect) %>%
  summarise(fle_stomach_content = sum(value)) %>% 
  arrange(name, year, ices_rect, quarter) %>%
  ungroup() %>% 
  mutate(id = paste(year, quarter, ices_rect, sep = "_")) %>% 
  group_by(id) %>%
  mutate(fle_stomach_content_tot = sum(fle_stomach_content),
         fle_stomach_content_prop = fle_stomach_content / fle_stomach_content_tot) %>% 
  ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 2180x32, now 34880x18]
#> group_by: 4 grouping variables (name, year, quarter, ices_rect)
#> summarise: now 1,296 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 81 unique values and 0% NA
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'fle_stomach_content_tot' (double) with 81 unique values and 0% NA
#>                   new variable 'fle_stomach_content_prop' (double) with 445 unique values and 1% NA
#> ungroup: no grouping variables

# This should amount to the number of unique prey
number_of_prey <- fle_prey_long %>% group_by(id) %>% mutate(n = n()) %>% ungroup() %>% distinct(n)
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'n' (integer) with one unique value and 0% NA
#> ungroup: no grouping variables
#> distinct: removed 1,295 rows (>99%), one row remaining
number_of_prey
#> # A tibble: 1 × 1
#>       n
#>   <int>
#> 1    16

test_id <- head(fle_prey_long$id, 1)
fle_prey_long %>% filter(id == test_id) %>% as.data.frame()
#> filter: removed 1,280 rows (99%), 16 rows remaining
#>                      name year quarter ices_rect fle_stomach_content
#> 1           amphipoda_tot 2015       4      40G4                1.30
#> 2            bivalvia_tot 2015       4      40G4              240.98
#> 3     clupea_harengus_tot 2015       4      40G4                0.00
#> 4           clupeidae_tot 2015       4      40G4                0.00
#> 5          gadiformes_tot 2015       4      40G4                0.00
#> 6        gadus_morhua_tot 2015       4      40G4                0.00
#> 7            gobiidae_tot 2015       4      40G4                0.00
#> 8             mysidae_tot 2015       4      40G4                0.00
#> 9             non_bio_tot 2015       4      40G4                5.78
#> 10    other_crustacea_tot 2015       4      40G4                2.12
#> 11       other_pisces_tot 2015       4      40G4                0.00
#> 12              other_tot 2015       4      40G4                2.27
#> 13 platichthys_flesus_tot 2015       4      40G4                0.00
#> 14         polychaeta_tot 2015       4      40G4                2.82
#> 15    saduria_entomon_tot 2015       4      40G4                0.07
#> 16  sprattus_sprattus_tot 2015       4      40G4                0.00
#>             id fle_stomach_content_tot fle_stomach_content_prop
#> 1  2015_4_40G4                  255.34             0.0050912509
#> 2  2015_4_40G4                  255.34             0.9437612595
#> 3  2015_4_40G4                  255.34             0.0000000000
#> 4  2015_4_40G4                  255.34             0.0000000000
#> 5  2015_4_40G4                  255.34             0.0000000000
#> 6  2015_4_40G4                  255.34             0.0000000000
#> 7  2015_4_40G4                  255.34             0.0000000000
#> 8  2015_4_40G4                  255.34             0.0000000000
#> 9  2015_4_40G4                  255.34             0.0226364847
#> 10 2015_4_40G4                  255.34             0.0083026553
#> 11 2015_4_40G4                  255.34             0.0000000000
#> 12 2015_4_40G4                  255.34             0.0088901073
#> 13 2015_4_40G4                  255.34             0.0000000000
#> 14 2015_4_40G4                  255.34             0.0110440981
#> 15 2015_4_40G4                  255.34             0.0002741443
#> 16 2015_4_40G4                  255.34             0.0000000000

# Now cod
# colnames(cod_prey)
cod_prey_long <- cod_prey %>%
  pivot_longer(15:30) %>% 
  group_by(name, year, quarter, ices_rect) %>%
  summarise(cod_stomach_content = sum(value)) %>% 
  arrange(name, year, ices_rect, quarter) %>%
  ungroup() %>% 
  mutate(id = paste(year, quarter, ices_rect, sep = "_")) %>% 
  group_by(id) %>%
  mutate(cod_stomach_content_tot = sum(cod_stomach_content),
         cod_stomach_content_prop = cod_stomach_content / cod_stomach_content_tot) %>% 
  ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x32, now 57488x18]
#> group_by: 4 grouping variables (name, year, quarter, ices_rect)
#> summarise: now 1,696 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 106 unique values and 0% NA
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'cod_stomach_content_tot' (double) with 106 unique values and 0% NA
#>                   new variable 'cod_stomach_content_prop' (double) with 784 unique values and 0% NA
#> ungroup: no grouping variables

unique(is.na(fle_prey_long))
#>       name  year quarter ices_rect fle_stomach_content    id
#> [1,] FALSE FALSE   FALSE     FALSE               FALSE FALSE
#> [2,] FALSE FALSE   FALSE     FALSE               FALSE FALSE
#>      fle_stomach_content_tot fle_stomach_content_prop
#> [1,]                   FALSE                    FALSE
#> [2,]                   FALSE                     TRUE
unique(is.na(cod_prey_long))
#>       name  year quarter ices_rect cod_stomach_content    id
#> [1,] FALSE FALSE   FALSE     FALSE               FALSE FALSE
#>      cod_stomach_content_tot cod_stomach_content_prop
#> [1,]                   FALSE                    FALSE

fle_prey_long %>% mutate(fle_stomach_content_prop = replace_na(fle_stomach_content_prop, -9)) %>% filter(fle_stomach_content_prop == -9)
#> mutate: changed 16 values (1%) of 'fle_stomach_content_prop' (16 fewer NA)
#> filter: removed 1,280 rows (99%), 16 rows remaining
#> # A tibble: 16 × 8
#>    name        year quarter ices_rect fle_stomach_cont… id     fle_stomach_cont…
#>    <chr>      <dbl>   <dbl> <chr>                 <dbl> <chr>              <dbl>
#>  1 amphipoda…  2016       1 40G8                      0 2016_…                 0
#>  2 bivalvia_…  2016       1 40G8                      0 2016_…                 0
#>  3 clupea_ha…  2016       1 40G8                      0 2016_…                 0
#>  4 clupeidae…  2016       1 40G8                      0 2016_…                 0
#>  5 gadiforme…  2016       1 40G8                      0 2016_…                 0
#>  6 gadus_mor…  2016       1 40G8                      0 2016_…                 0
#>  7 gobiidae_…  2016       1 40G8                      0 2016_…                 0
#>  8 mysidae_t…  2016       1 40G8                      0 2016_…                 0
#>  9 non_bio_t…  2016       1 40G8                      0 2016_…                 0
#> 10 other_cru…  2016       1 40G8                      0 2016_…                 0
#> 11 other_pis…  2016       1 40G8                      0 2016_…                 0
#> 12 other_tot   2016       1 40G8                      0 2016_…                 0
#> 13 platichth…  2016       1 40G8                      0 2016_…                 0
#> 14 polychaet…  2016       1 40G8                      0 2016_…                 0
#> 15 saduria_e…  2016       1 40G8                      0 2016_…                 0
#> 16 sprattus_…  2016       1 40G8                      0 2016_…                 0
#> # … with 1 more variable: fle_stomach_content_prop <dbl>
fle_prey_long <- fle_prey_long %>% mutate(fle_stomach_content_prop = replace_na(fle_stomach_content_prop, 0)) 
#> mutate: changed 16 values (1%) of 'fle_stomach_content_prop' (16 fewer NA)

# Calculate Schoener index
schoener <- left_join(cod_prey_long, fle_prey_long) %>% 
  drop_na(name) %>% 
  drop_na(fle_stomach_content_prop) %>% 
  drop_na(cod_stomach_content_prop) %>% 
  group_by(year, quarter, ices_rect) %>%
  summarise(schoener = 1 - 0.5*(sum(abs(fle_stomach_content_prop - cod_stomach_content_prop)))) %>% 
  ungroup()
#> Joining, by = c("name", "year", "quarter", "ices_rect", "id")
#> left_join: added 3 columns (fle_stomach_content, fle_stomach_content_tot, fle_stomach_content_prop)
#>            > rows only in x     432
#>            > rows only in y  (   32)
#>            > matched rows     1,264
#>            >                 =======
#>            > rows total       1,696
#> drop_na: no rows removed
#> drop_na: removed 432 rows (25%), 1,264 rows remaining
#> drop_na: no rows removed
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> summarise: now 79 rows and 4 columns, 2 group variables remaining (year, quarter)
#> ungroup: no grouping variables

# Summarise cod and flounder data by ices_rect then add to diet data
colnames(cod)
#>  [1] "FR_tot"                   "FR_sad"                  
#>  [3] "FR_spr"                   "FR_her"                  
#>  [5] "saduria_entomon_tot"      "tot_prey_biom"           
#>  [7] "common_tot"               "unique_pred_id"          
#>  [9] "year"                     "quarter"                 
#> [11] "time_period"              "quarter_fact"            
#> [13] "pred_weight_g"            "pred_weight_source"      
#> [15] "pred_cm"                  "predator_spec"           
#> [17] "predcod_density"          "predfle_density"         
#> [19] "predcod_density_sc"       "predfle_density_sc"      
#> [21] "depth"                    "X"                       
#> [23] "Y"                        "lat"                     
#> [25] "long"                     "ices_rect"               
#> [27] "sub_div"                  "cruise"                  
#> [29] "depth2_sc"                "saduria_entomon_per_mass"
#> [31] "tot_prey_biom_per_mass"   "depth_sc"
cod$year_rect_id <- paste(cod$year, cod$quarter, cod$ices_rect, sep = "_")

dens_sum <- cod %>%
  drop_na(predfle_density) %>% 
  drop_na(predcod_density) %>% 
  group_by(year_rect_id) %>%
  summarise(predfle_density_tot = sum(predfle_density),
            predcod_density_tot = sum(predcod_density)) %>% 
  ungroup() %>% 
  mutate(predfle_density_tot_sc = (predfle_density_tot - mean(predfle_density_tot)) / sd(predfle_density_tot),
         predcod_density_tot_sc = (predcod_density_tot - mean(predcod_density_tot)) / sd(predcod_density_tot))
#> drop_na: no rows removed
#> drop_na: no rows removed
#> group_by: one grouping variable (year_rect_id)
#> summarise: now 106 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'predfle_density_tot_sc' (double) with 106 unique values and 0% NA
#>         new variable 'predcod_density_tot_sc' (double) with 106 unique values and 0% NA

schoener$year_rect_id <- paste(schoener$year, schoener$quarter, schoener$ices_rect, sep = "_")

schoener <- left_join(schoener, dens_sum)
#> Joining, by = "year_rect_id"
#> left_join: added 4 columns (predfle_density_tot, predcod_density_tot, predfle_density_tot_sc, predcod_density_tot_sc)
#>            > rows only in x    0
#>            > rows only in y  (27)
#>            > matched rows     79
#>            >                 ====
#>            > rows total       79

# Quickly check data to determine which distribution to use
schoener %>% 
  ungroup() %>% 
  count(schoener == 0) %>% 
  mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now 2 rows and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with 2 unique values and 0% NA
#> # A tibble: 2 × 3
#>   `schoener == 0`     n   prop
#>   <lgl>           <int>  <dbl>
#> 1 FALSE              76 0.962 
#> 2 TRUE                3 0.0380

# Fit beta models, so few zeroes and no 1's
schoener %>% arrange(desc(schoener)) %>% dplyr::select(schoener)
#> # A tibble: 79 × 1
#>    schoener
#>       <dbl>
#>  1    0.624
#>  2    0.614
#>  3    0.565
#>  4    0.556
#>  5    0.553
#>  6    0.5  
#>  7    0.349
#>  8    0.308
#>  9    0.222
#> 10    0.208
#> # … with 69 more rows

# Final polish of data before feeding into models (species, not size-based indicies)
schoener <- schoener %>%
  mutate(schoener2 = ifelse(schoener == 0, 0.0001, schoener),
         year_f = as.factor(year),
         quarter_f = as.factor(quarter),
         ices_rect_f = as.factor(ices_rect))
#> mutate: new variable 'schoener2' (double) with 77 unique values and 0% NA
#>         new variable 'year_f' (factor) with 6 unique values and 0% NA
#>         new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#>         new variable 'ices_rect_f' (factor) with 18 unique values and 0% NA

Calculate Schoener’s overlap index for diet per size groups (cod 0-30, 30+ and all flounder sizes)

# colnames(cod_prey)
cod_prey_long2 <- cod_prey %>%
  mutate(size_group = ifelse(pred_cm > 30, "Large", "Small")) %>% 
  pivot_longer(15:30) %>% 
  group_by(name, year, quarter, ices_rect, size_group) %>%
  summarise(cod_stomach_content = sum(value)) %>% 
  arrange(name, year, ices_rect, quarter) %>%
  ungroup() %>% 
  mutate(id = paste(year, quarter, ices_rect, sep = "_")) %>% 
  group_by(id, size_group) %>% # Add size-group here
  mutate(cod_stomach_content_tot = sum(cod_stomach_content),
         cod_stomach_content_prop = cod_stomach_content
         / cod_stomach_content_tot) %>% 
  ungroup()
#> mutate: new variable 'size_group' (character) with 2 unique values and 0% NA
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x33, now 57488x19]
#> group_by: 5 grouping variables (name, year, quarter, ices_rect, size_group)
#> summarise: now 3,200 rows and 6 columns, 4 group variables remaining (name, year, quarter, ices_rect)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 106 unique values and 0% NA
#> group_by: 2 grouping variables (id, size_group)
#> mutate (grouped): new variable 'cod_stomach_content_tot' (double) with 193 unique values and 0% NA
#>                   new variable 'cod_stomach_content_prop' (double) with 1,079 unique values and <1% NA
#> ungroup: no grouping variables

# This should amount to the number of unique prey
cod_prey_long2 %>% group_by(id, size_group) %>% mutate(n = n()) %>% ungroup() %>% distinct(n)
#> group_by: 2 grouping variables (id, size_group)
#> mutate (grouped): new variable 'n' (integer) with one unique value and 0% NA
#> ungroup: no grouping variables
#> distinct: removed 3,199 rows (>99%), one row remaining
#> # A tibble: 1 × 1
#>       n
#>   <int>
#> 1    16

# Split by size-group
cod_prey_long_small <- cod_prey_long2 %>% filter(size_group == "Small") %>%
  rename("cod_stomach_content_prop_small" = "cod_stomach_content_prop") %>%
  dplyr::select(name, id, cod_stomach_content_prop_small)
#> filter: removed 1,632 rows (51%), 1,568 rows remaining
#> rename: renamed one variable (cod_stomach_content_prop_small)

cod_prey_long_large <- cod_prey_long2 %>% filter(size_group == "Large") %>% 
  rename("cod_stomach_content_prop_large" = "cod_stomach_content_prop") %>%
  dplyr::select(name, id, cod_stomach_content_prop_large)
#> filter: removed 1,568 rows (49%), 1,632 rows remaining
#> rename: renamed one variable (cod_stomach_content_prop_large)

cod_prey_long_large %>% filter(cod_stomach_content_prop_large == "NaN")
#> filter: removed all rows (100%)
#> # A tibble: 0 × 3
#> # … with 3 variables: name <chr>, id <chr>,
#> #   cod_stomach_content_prop_large <dbl>
cod_prey_long_small %>% filter(cod_stomach_content_prop_small == "NaN")
#> filter: removed 1,552 rows (99%), 16 rows remaining
#> # A tibble: 16 × 3
#>    name                   id          cod_stomach_content_prop_small
#>    <chr>                  <chr>                                <dbl>
#>  1 amphipoda_tot          2017_1_43G9                            NaN
#>  2 bivalvia_tot           2017_1_43G9                            NaN
#>  3 clupea_harengus_tot    2017_1_43G9                            NaN
#>  4 clupeidae_tot          2017_1_43G9                            NaN
#>  5 gadiformes_tot         2017_1_43G9                            NaN
#>  6 gadus_morhua_tot       2017_1_43G9                            NaN
#>  7 gobiidae_tot           2017_1_43G9                            NaN
#>  8 mysidae_tot            2017_1_43G9                            NaN
#>  9 non_bio_tot            2017_1_43G9                            NaN
#> 10 other_crustacea_tot    2017_1_43G9                            NaN
#> 11 other_pisces_tot       2017_1_43G9                            NaN
#> 12 other_tot              2017_1_43G9                            NaN
#> 13 platichthys_flesus_tot 2017_1_43G9                            NaN
#> 14 polychaeta_tot         2017_1_43G9                            NaN
#> 15 saduria_entomon_tot    2017_1_43G9                            NaN
#> 16 sprattus_sprattus_tot  2017_1_43G9                            NaN
cod_prey_long_small <- cod_prey_long_small %>%
  mutate(cod_stomach_content_prop_small = ifelse(cod_stomach_content_prop_small == "NaN",
                                                 0,
                                                 cod_stomach_content_prop_small))
#> mutate: changed 16 values (1%) of 'cod_stomach_content_prop_small' (16 fewer NA)

# Calculate Schoener index
schoener2 <- fle_prey_long %>%
  left_join(cod_prey_long_small) %>% 
  left_join(cod_prey_long_large) %>% 
  drop_na(name) %>% 
  drop_na(fle_stomach_content_prop) %>% 
  drop_na(cod_stomach_content_prop_small) %>%
  drop_na(cod_stomach_content_prop_large) %>% 
  group_by(year, quarter, ices_rect) %>%
  summarise(schoener_f_sc = 1 - 0.5*(sum(abs(fle_stomach_content_prop - cod_stomach_content_prop_small))),
            schoener_f_lc = 1 - 0.5*(sum(abs(fle_stomach_content_prop - cod_stomach_content_prop_large))),
            schoener_sc_lc = 1 - 0.5*(sum(abs(cod_stomach_content_prop_small - cod_stomach_content_prop_large)))) %>% 
  ungroup()
#> Joining, by = c("name", "id")
#> left_join: added one column (cod_stomach_content_prop_small)
#>            > rows only in x      64
#>            > rows only in y  (  336)
#>            > matched rows     1,232
#>            >                 =======
#>            > rows total       1,296
#> Joining, by = c("name", "id")
#> left_join: added one column (cod_stomach_content_prop_large)
#>            > rows only in x      48
#>            > rows only in y  (  384)
#>            > matched rows     1,248
#>            >                 =======
#>            > rows total       1,296
#> drop_na: no rows removed
#> drop_na: no rows removed
#> drop_na: removed 64 rows (5%), 1,232 rows remaining
#> drop_na: removed 16 rows (1%), 1,216 rows remaining
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> summarise: now 76 rows and 6 columns, 2 group variables remaining (year, quarter)
#> ungroup: no grouping variables

schoener_long2 <- schoener2 %>%
  pivot_longer(4:6, names_to = "schoener_combination", values_to = "schoener")
#> pivot_longer: reorganized (schoener_f_sc, schoener_f_lc, schoener_sc_lc) into (schoener_combination, schoener) [was 76x6, now 228x5]

ggplot(schoener_long2, aes(schoener_combination, schoener, fill = factor(schoener_combination),
                           color = factor(schoener_combination))) +
  ggdist::stat_halfeye(adjust = 0.5, justification = -0.1, .width = 0, point_colour = NA, alpha = 0.8,
                       show.legend = FALSE) +
  geom_boxplot(width = 0.12, outlier.color = NA, alpha = 0.5, show.legend = FALSE) +
  ggdist::stat_dots(side = "left", justification = 1.1, alpha = 0.6) + 
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "",
                    labels = c("Flounder-Cod (L)", "Flounder-Cod (S)", "Cod (S)-Cod (L)")) + 
  guides(color = FALSE, fill = guide_legend(override.aes = list(
    shape = 21, size = 2, fill = brewer.pal(n = 3, name = "Dark2"), color = brewer.pal(n = 3, name = "Dark2")))) +
  coord_flip() + 
  labs(y = "Value", x = "") +
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        axis.text.y = element_blank()) +
  ggtitle("Schoeners's overlap index") + 
  NULL
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.

ggsave("figures/schoener_size_groups_data.png", width = 6.5, height = 6.5, dpi = 600)

# Prepare for analysis
schoener_long2 <- schoener_long2 %>%
  mutate(schoener2 = ifelse(schoener == 0, 0.0001, schoener),
         schoener2 = ifelse(schoener == 1, 1-0.0001, schoener2),
         year_f = as.factor(year),
         quarter_f = as.factor(quarter),
         ices_rect_f = as.factor(ices_rect),
         schoener_combination_f = as.factor(schoener_combination))
#> mutate: new variable 'schoener2' (double) with 204 unique values and 0% NA
#>         new variable 'year_f' (factor) with 6 unique values and 0% NA
#>         new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#>         new variable 'ices_rect_f' (factor) with 18 unique values and 0% NA
#>         new variable 'schoener_combination_f' (factor) with 3 unique values and 0% NA

Fit brms models of diversity and overlap indices

Beta model schoener index with density covariates

# fit
m_schoen_beta <- brm(
  bf(schoener2 ~ 0 + year_f + quarter_f + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
     phi ~  0 + year_f + quarter_f + predfle_density_tot_sc + predcod_density_tot_sc),
  data = schoener, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
  chains = 4, iter = 4000, cores = 4, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling

plot(m_schoen_beta)

loo_m_schoen_beta <- loo(m_schoen_beta, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.

Plot

# Evaluate fit and convergence etc.
# PP check
pp_check(m_schoen_beta, ndraws = 50) +
  theme_light(base_size = 20) + 
  scale_color_brewer(palette = "Dark2", name = "") +
  NULL
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.

ggsave("figures/supp/schoener_pp_check.png", width = 6.5, height = 6.5, dpi = 600)

# Chain convergence
posterior <- as.array(m_schoen_beta)
dimnames(posterior)
#> $iteration
#>    [1] "1"    "2"    "3"    "4"    "5"    "6"    "7"    "8"    "9"    "10"  
#>   [11] "11"   "12"   "13"   "14"   "15"   "16"   "17"   "18"   "19"   "20"  
#>   [21] "21"   "22"   "23"   "24"   "25"   "26"   "27"   "28"   "29"   "30"  
#>   [31] "31"   "32"   "33"   "34"   "35"   "36"   "37"   "38"   "39"   "40"  
#>   [41] "41"   "42"   "43"   "44"   "45"   "46"   "47"   "48"   "49"   "50"  
#>   [51] "51"   "52"   "53"   "54"   "55"   "56"   "57"   "58"   "59"   "60"  
#>   [61] "61"   "62"   "63"   "64"   "65"   "66"   "67"   "68"   "69"   "70"  
#>   [71] "71"   "72"   "73"   "74"   "75"   "76"   "77"   "78"   "79"   "80"  
#>   [81] "81"   "82"   "83"   "84"   "85"   "86"   "87"   "88"   "89"   "90"  
#>   [91] "91"   "92"   "93"   "94"   "95"   "96"   "97"   "98"   "99"   "100" 
#>  [101] "101"  "102"  "103"  "104"  "105"  "106"  "107"  "108"  "109"  "110" 
#>  [111] "111"  "112"  "113"  "114"  "115"  "116"  "117"  "118"  "119"  "120" 
#>  [121] "121"  "122"  "123"  "124"  "125"  "126"  "127"  "128"  "129"  "130" 
#>  [131] "131"  "132"  "133"  "134"  "135"  "136"  "137"  "138"  "139"  "140" 
#>  [141] "141"  "142"  "143"  "144"  "145"  "146"  "147"  "148"  "149"  "150" 
#>  [151] "151"  "152"  "153"  "154"  "155"  "156"  "157"  "158"  "159"  "160" 
#>  [161] "161"  "162"  "163"  "164"  "165"  "166"  "167"  "168"  "169"  "170" 
#>  [171] "171"  "172"  "173"  "174"  "175"  "176"  "177"  "178"  "179"  "180" 
#>  [181] "181"  "182"  "183"  "184"  "185"  "186"  "187"  "188"  "189"  "190" 
#>  [191] "191"  "192"  "193"  "194"  "195"  "196"  "197"  "198"  "199"  "200" 
#>  [201] "201"  "202"  "203"  "204"  "205"  "206"  "207"  "208"  "209"  "210" 
#>  [211] "211"  "212"  "213"  "214"  "215"  "216"  "217"  "218"  "219"  "220" 
#>  [221] "221"  "222"  "223"  "224"  "225"  "226"  "227"  "228"  "229"  "230" 
#>  [231] "231"  "232"  "233"  "234"  "235"  "236"  "237"  "238"  "239"  "240" 
#>  [241] "241"  "242"  "243"  "244"  "245"  "246"  "247"  "248"  "249"  "250" 
#>  [251] "251"  "252"  "253"  "254"  "255"  "256"  "257"  "258"  "259"  "260" 
#>  [261] "261"  "262"  "263"  "264"  "265"  "266"  "267"  "268"  "269"  "270" 
#>  [271] "271"  "272"  "273"  "274"  "275"  "276"  "277"  "278"  "279"  "280" 
#>  [281] "281"  "282"  "283"  "284"  "285"  "286"  "287"  "288"  "289"  "290" 
#>  [291] "291"  "292"  "293"  "294"  "295"  "296"  "297"  "298"  "299"  "300" 
#>  [301] "301"  "302"  "303"  "304"  "305"  "306"  "307"  "308"  "309"  "310" 
#>  [311] "311"  "312"  "313"  "314"  "315"  "316"  "317"  "318"  "319"  "320" 
#>  [321] "321"  "322"  "323"  "324"  "325"  "326"  "327"  "328"  "329"  "330" 
#>  [331] "331"  "332"  "333"  "334"  "335"  "336"  "337"  "338"  "339"  "340" 
#>  [341] "341"  "342"  "343"  "344"  "345"  "346"  "347"  "348"  "349"  "350" 
#>  [351] "351"  "352"  "353"  "354"  "355"  "356"  "357"  "358"  "359"  "360" 
#>  [361] "361"  "362"  "363"  "364"  "365"  "366"  "367"  "368"  "369"  "370" 
#>  [371] "371"  "372"  "373"  "374"  "375"  "376"  "377"  "378"  "379"  "380" 
#>  [381] "381"  "382"  "383"  "384"  "385"  "386"  "387"  "388"  "389"  "390" 
#>  [391] "391"  "392"  "393"  "394"  "395"  "396"  "397"  "398"  "399"  "400" 
#>  [401] "401"  "402"  "403"  "404"  "405"  "406"  "407"  "408"  "409"  "410" 
#>  [411] "411"  "412"  "413"  "414"  "415"  "416"  "417"  "418"  "419"  "420" 
#>  [421] "421"  "422"  "423"  "424"  "425"  "426"  "427"  "428"  "429"  "430" 
#>  [431] "431"  "432"  "433"  "434"  "435"  "436"  "437"  "438"  "439"  "440" 
#>  [441] "441"  "442"  "443"  "444"  "445"  "446"  "447"  "448"  "449"  "450" 
#>  [451] "451"  "452"  "453"  "454"  "455"  "456"  "457"  "458"  "459"  "460" 
#>  [461] "461"  "462"  "463"  "464"  "465"  "466"  "467"  "468"  "469"  "470" 
#>  [471] "471"  "472"  "473"  "474"  "475"  "476"  "477"  "478"  "479"  "480" 
#>  [481] "481"  "482"  "483"  "484"  "485"  "486"  "487"  "488"  "489"  "490" 
#>  [491] "491"  "492"  "493"  "494"  "495"  "496"  "497"  "498"  "499"  "500" 
#>  [501] "501"  "502"  "503"  "504"  "505"  "506"  "507"  "508"  "509"  "510" 
#>  [511] "511"  "512"  "513"  "514"  "515"  "516"  "517"  "518"  "519"  "520" 
#>  [521] "521"  "522"  "523"  "524"  "525"  "526"  "527"  "528"  "529"  "530" 
#>  [531] "531"  "532"  "533"  "534"  "535"  "536"  "537"  "538"  "539"  "540" 
#>  [541] "541"  "542"  "543"  "544"  "545"  "546"  "547"  "548"  "549"  "550" 
#>  [551] "551"  "552"  "553"  "554"  "555"  "556"  "557"  "558"  "559"  "560" 
#>  [561] "561"  "562"  "563"  "564"  "565"  "566"  "567"  "568"  "569"  "570" 
#>  [571] "571"  "572"  "573"  "574"  "575"  "576"  "577"  "578"  "579"  "580" 
#>  [581] "581"  "582"  "583"  "584"  "585"  "586"  "587"  "588"  "589"  "590" 
#>  [591] "591"  "592"  "593"  "594"  "595"  "596"  "597"  "598"  "599"  "600" 
#>  [601] "601"  "602"  "603"  "604"  "605"  "606"  "607"  "608"  "609"  "610" 
#>  [611] "611"  "612"  "613"  "614"  "615"  "616"  "617"  "618"  "619"  "620" 
#>  [621] "621"  "622"  "623"  "624"  "625"  "626"  "627"  "628"  "629"  "630" 
#>  [631] "631"  "632"  "633"  "634"  "635"  "636"  "637"  "638"  "639"  "640" 
#>  [641] "641"  "642"  "643"  "644"  "645"  "646"  "647"  "648"  "649"  "650" 
#>  [651] "651"  "652"  "653"  "654"  "655"  "656"  "657"  "658"  "659"  "660" 
#>  [661] "661"  "662"  "663"  "664"  "665"  "666"  "667"  "668"  "669"  "670" 
#>  [671] "671"  "672"  "673"  "674"  "675"  "676"  "677"  "678"  "679"  "680" 
#>  [681] "681"  "682"  "683"  "684"  "685"  "686"  "687"  "688"  "689"  "690" 
#>  [691] "691"  "692"  "693"  "694"  "695"  "696"  "697"  "698"  "699"  "700" 
#>  [701] "701"  "702"  "703"  "704"  "705"  "706"  "707"  "708"  "709"  "710" 
#>  [711] "711"  "712"  "713"  "714"  "715"  "716"  "717"  "718"  "719"  "720" 
#>  [721] "721"  "722"  "723"  "724"  "725"  "726"  "727"  "728"  "729"  "730" 
#>  [731] "731"  "732"  "733"  "734"  "735"  "736"  "737"  "738"  "739"  "740" 
#>  [741] "741"  "742"  "743"  "744"  "745"  "746"  "747"  "748"  "749"  "750" 
#>  [751] "751"  "752"  "753"  "754"  "755"  "756"  "757"  "758"  "759"  "760" 
#>  [761] "761"  "762"  "763"  "764"  "765"  "766"  "767"  "768"  "769"  "770" 
#>  [771] "771"  "772"  "773"  "774"  "775"  "776"  "777"  "778"  "779"  "780" 
#>  [781] "781"  "782"  "783"  "784"  "785"  "786"  "787"  "788"  "789"  "790" 
#>  [791] "791"  "792"  "793"  "794"  "795"  "796"  "797"  "798"  "799"  "800" 
#>  [801] "801"  "802"  "803"  "804"  "805"  "806"  "807"  "808"  "809"  "810" 
#>  [811] "811"  "812"  "813"  "814"  "815"  "816"  "817"  "818"  "819"  "820" 
#>  [821] "821"  "822"  "823"  "824"  "825"  "826"  "827"  "828"  "829"  "830" 
#>  [831] "831"  "832"  "833"  "834"  "835"  "836"  "837"  "838"  "839"  "840" 
#>  [841] "841"  "842"  "843"  "844"  "845"  "846"  "847"  "848"  "849"  "850" 
#>  [851] "851"  "852"  "853"  "854"  "855"  "856"  "857"  "858"  "859"  "860" 
#>  [861] "861"  "862"  "863"  "864"  "865"  "866"  "867"  "868"  "869"  "870" 
#>  [871] "871"  "872"  "873"  "874"  "875"  "876"  "877"  "878"  "879"  "880" 
#>  [881] "881"  "882"  "883"  "884"  "885"  "886"  "887"  "888"  "889"  "890" 
#>  [891] "891"  "892"  "893"  "894"  "895"  "896"  "897"  "898"  "899"  "900" 
#>  [901] "901"  "902"  "903"  "904"  "905"  "906"  "907"  "908"  "909"  "910" 
#>  [911] "911"  "912"  "913"  "914"  "915"  "916"  "917"  "918"  "919"  "920" 
#>  [921] "921"  "922"  "923"  "924"  "925"  "926"  "927"  "928"  "929"  "930" 
#>  [931] "931"  "932"  "933"  "934"  "935"  "936"  "937"  "938"  "939"  "940" 
#>  [941] "941"  "942"  "943"  "944"  "945"  "946"  "947"  "948"  "949"  "950" 
#>  [951] "951"  "952"  "953"  "954"  "955"  "956"  "957"  "958"  "959"  "960" 
#>  [961] "961"  "962"  "963"  "964"  "965"  "966"  "967"  "968"  "969"  "970" 
#>  [971] "971"  "972"  "973"  "974"  "975"  "976"  "977"  "978"  "979"  "980" 
#>  [981] "981"  "982"  "983"  "984"  "985"  "986"  "987"  "988"  "989"  "990" 
#>  [991] "991"  "992"  "993"  "994"  "995"  "996"  "997"  "998"  "999"  "1000"
#> [1001] "1001" "1002" "1003" "1004" "1005" "1006" "1007" "1008" "1009" "1010"
#> [1011] "1011" "1012" "1013" "1014" "1015" "1016" "1017" "1018" "1019" "1020"
#> [1021] "1021" "1022" "1023" "1024" "1025" "1026" "1027" "1028" "1029" "1030"
#> [1031] "1031" "1032" "1033" "1034" "1035" "1036" "1037" "1038" "1039" "1040"
#> [1041] "1041" "1042" "1043" "1044" "1045" "1046" "1047" "1048" "1049" "1050"
#> [1051] "1051" "1052" "1053" "1054" "1055" "1056" "1057" "1058" "1059" "1060"
#> [1061] "1061" "1062" "1063" "1064" "1065" "1066" "1067" "1068" "1069" "1070"
#> [1071] "1071" "1072" "1073" "1074" "1075" "1076" "1077" "1078" "1079" "1080"
#> [1081] "1081" "1082" "1083" "1084" "1085" "1086" "1087" "1088" "1089" "1090"
#> [1091] "1091" "1092" "1093" "1094" "1095" "1096" "1097" "1098" "1099" "1100"
#> [1101] "1101" "1102" "1103" "1104" "1105" "1106" "1107" "1108" "1109" "1110"
#> [1111] "1111" "1112" "1113" "1114" "1115" "1116" "1117" "1118" "1119" "1120"
#> [1121] "1121" "1122" "1123" "1124" "1125" "1126" "1127" "1128" "1129" "1130"
#> [1131] "1131" "1132" "1133" "1134" "1135" "1136" "1137" "1138" "1139" "1140"
#> [1141] "1141" "1142" "1143" "1144" "1145" "1146" "1147" "1148" "1149" "1150"
#> [1151] "1151" "1152" "1153" "1154" "1155" "1156" "1157" "1158" "1159" "1160"
#> [1161] "1161" "1162" "1163" "1164" "1165" "1166" "1167" "1168" "1169" "1170"
#> [1171] "1171" "1172" "1173" "1174" "1175" "1176" "1177" "1178" "1179" "1180"
#> [1181] "1181" "1182" "1183" "1184" "1185" "1186" "1187" "1188" "1189" "1190"
#> [1191] "1191" "1192" "1193" "1194" "1195" "1196" "1197" "1198" "1199" "1200"
#> [1201] "1201" "1202" "1203" "1204" "1205" "1206" "1207" "1208" "1209" "1210"
#> [1211] "1211" "1212" "1213" "1214" "1215" "1216" "1217" "1218" "1219" "1220"
#> [1221] "1221" "1222" "1223" "1224" "1225" "1226" "1227" "1228" "1229" "1230"
#> [1231] "1231" "1232" "1233" "1234" "1235" "1236" "1237" "1238" "1239" "1240"
#> [1241] "1241" "1242" "1243" "1244" "1245" "1246" "1247" "1248" "1249" "1250"
#> [1251] "1251" "1252" "1253" "1254" "1255" "1256" "1257" "1258" "1259" "1260"
#> [1261] "1261" "1262" "1263" "1264" "1265" "1266" "1267" "1268" "1269" "1270"
#> [1271] "1271" "1272" "1273" "1274" "1275" "1276" "1277" "1278" "1279" "1280"
#> [1281] "1281" "1282" "1283" "1284" "1285" "1286" "1287" "1288" "1289" "1290"
#> [1291] "1291" "1292" "1293" "1294" "1295" "1296" "1297" "1298" "1299" "1300"
#> [1301] "1301" "1302" "1303" "1304" "1305" "1306" "1307" "1308" "1309" "1310"
#> [1311] "1311" "1312" "1313" "1314" "1315" "1316" "1317" "1318" "1319" "1320"
#> [1321] "1321" "1322" "1323" "1324" "1325" "1326" "1327" "1328" "1329" "1330"
#> [1331] "1331" "1332" "1333" "1334" "1335" "1336" "1337" "1338" "1339" "1340"
#> [1341] "1341" "1342" "1343" "1344" "1345" "1346" "1347" "1348" "1349" "1350"
#> [1351] "1351" "1352" "1353" "1354" "1355" "1356" "1357" "1358" "1359" "1360"
#> [1361] "1361" "1362" "1363" "1364" "1365" "1366" "1367" "1368" "1369" "1370"
#> [1371] "1371" "1372" "1373" "1374" "1375" "1376" "1377" "1378" "1379" "1380"
#> [1381] "1381" "1382" "1383" "1384" "1385" "1386" "1387" "1388" "1389" "1390"
#> [1391] "1391" "1392" "1393" "1394" "1395" "1396" "1397" "1398" "1399" "1400"
#> [1401] "1401" "1402" "1403" "1404" "1405" "1406" "1407" "1408" "1409" "1410"
#> [1411] "1411" "1412" "1413" "1414" "1415" "1416" "1417" "1418" "1419" "1420"
#> [1421] "1421" "1422" "1423" "1424" "1425" "1426" "1427" "1428" "1429" "1430"
#> [1431] "1431" "1432" "1433" "1434" "1435" "1436" "1437" "1438" "1439" "1440"
#> [1441] "1441" "1442" "1443" "1444" "1445" "1446" "1447" "1448" "1449" "1450"
#> [1451] "1451" "1452" "1453" "1454" "1455" "1456" "1457" "1458" "1459" "1460"
#> [1461] "1461" "1462" "1463" "1464" "1465" "1466" "1467" "1468" "1469" "1470"
#> [1471] "1471" "1472" "1473" "1474" "1475" "1476" "1477" "1478" "1479" "1480"
#> [1481] "1481" "1482" "1483" "1484" "1485" "1486" "1487" "1488" "1489" "1490"
#> [1491] "1491" "1492" "1493" "1494" "1495" "1496" "1497" "1498" "1499" "1500"
#> [1501] "1501" "1502" "1503" "1504" "1505" "1506" "1507" "1508" "1509" "1510"
#> [1511] "1511" "1512" "1513" "1514" "1515" "1516" "1517" "1518" "1519" "1520"
#> [1521] "1521" "1522" "1523" "1524" "1525" "1526" "1527" "1528" "1529" "1530"
#> [1531] "1531" "1532" "1533" "1534" "1535" "1536" "1537" "1538" "1539" "1540"
#> [1541] "1541" "1542" "1543" "1544" "1545" "1546" "1547" "1548" "1549" "1550"
#> [1551] "1551" "1552" "1553" "1554" "1555" "1556" "1557" "1558" "1559" "1560"
#> [1561] "1561" "1562" "1563" "1564" "1565" "1566" "1567" "1568" "1569" "1570"
#> [1571] "1571" "1572" "1573" "1574" "1575" "1576" "1577" "1578" "1579" "1580"
#> [1581] "1581" "1582" "1583" "1584" "1585" "1586" "1587" "1588" "1589" "1590"
#> [1591] "1591" "1592" "1593" "1594" "1595" "1596" "1597" "1598" "1599" "1600"
#> [1601] "1601" "1602" "1603" "1604" "1605" "1606" "1607" "1608" "1609" "1610"
#> [1611] "1611" "1612" "1613" "1614" "1615" "1616" "1617" "1618" "1619" "1620"
#> [1621] "1621" "1622" "1623" "1624" "1625" "1626" "1627" "1628" "1629" "1630"
#> [1631] "1631" "1632" "1633" "1634" "1635" "1636" "1637" "1638" "1639" "1640"
#> [1641] "1641" "1642" "1643" "1644" "1645" "1646" "1647" "1648" "1649" "1650"
#> [1651] "1651" "1652" "1653" "1654" "1655" "1656" "1657" "1658" "1659" "1660"
#> [1661] "1661" "1662" "1663" "1664" "1665" "1666" "1667" "1668" "1669" "1670"
#> [1671] "1671" "1672" "1673" "1674" "1675" "1676" "1677" "1678" "1679" "1680"
#> [1681] "1681" "1682" "1683" "1684" "1685" "1686" "1687" "1688" "1689" "1690"
#> [1691] "1691" "1692" "1693" "1694" "1695" "1696" "1697" "1698" "1699" "1700"
#> [1701] "1701" "1702" "1703" "1704" "1705" "1706" "1707" "1708" "1709" "1710"
#> [1711] "1711" "1712" "1713" "1714" "1715" "1716" "1717" "1718" "1719" "1720"
#> [1721] "1721" "1722" "1723" "1724" "1725" "1726" "1727" "1728" "1729" "1730"
#> [1731] "1731" "1732" "1733" "1734" "1735" "1736" "1737" "1738" "1739" "1740"
#> [1741] "1741" "1742" "1743" "1744" "1745" "1746" "1747" "1748" "1749" "1750"
#> [1751] "1751" "1752" "1753" "1754" "1755" "1756" "1757" "1758" "1759" "1760"
#> [1761] "1761" "1762" "1763" "1764" "1765" "1766" "1767" "1768" "1769" "1770"
#> [1771] "1771" "1772" "1773" "1774" "1775" "1776" "1777" "1778" "1779" "1780"
#> [1781] "1781" "1782" "1783" "1784" "1785" "1786" "1787" "1788" "1789" "1790"
#> [1791] "1791" "1792" "1793" "1794" "1795" "1796" "1797" "1798" "1799" "1800"
#> [1801] "1801" "1802" "1803" "1804" "1805" "1806" "1807" "1808" "1809" "1810"
#> [1811] "1811" "1812" "1813" "1814" "1815" "1816" "1817" "1818" "1819" "1820"
#> [1821] "1821" "1822" "1823" "1824" "1825" "1826" "1827" "1828" "1829" "1830"
#> [1831] "1831" "1832" "1833" "1834" "1835" "1836" "1837" "1838" "1839" "1840"
#> [1841] "1841" "1842" "1843" "1844" "1845" "1846" "1847" "1848" "1849" "1850"
#> [1851] "1851" "1852" "1853" "1854" "1855" "1856" "1857" "1858" "1859" "1860"
#> [1861] "1861" "1862" "1863" "1864" "1865" "1866" "1867" "1868" "1869" "1870"
#> [1871] "1871" "1872" "1873" "1874" "1875" "1876" "1877" "1878" "1879" "1880"
#> [1881] "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889" "1890"
#> [1891] "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899" "1900"
#> [1901] "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909" "1910"
#> [1911] "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919" "1920"
#> [1921] "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929" "1930"
#> [1931] "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939" "1940"
#> [1941] "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949" "1950"
#> [1951] "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959" "1960"
#> [1961] "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969" "1970"
#> [1971] "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979" "1980"
#> [1981] "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989" "1990"
#> [1991] "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999" "2000"
#> 
#> $chain
#> [1] "1" "2" "3" "4"
#> 
#> $variable
#>  [1] "b_year_f2015"                  "b_year_f2016"                 
#>  [3] "b_year_f2017"                  "b_year_f2018"                 
#>  [5] "b_year_f2019"                  "b_year_f2020"                 
#>  [7] "b_quarter_f4"                  "b_predfle_density_tot_sc"     
#>  [9] "b_predcod_density_tot_sc"      "b_phi_year_f2015"             
#> [11] "b_phi_year_f2016"              "b_phi_year_f2017"             
#> [13] "b_phi_year_f2018"              "b_phi_year_f2019"             
#> [15] "b_phi_year_f2020"              "b_phi_quarter_f4"             
#> [17] "b_phi_predfle_density_tot_sc"  "b_phi_predcod_density_tot_sc" 
#> [19] "sd_ices_rect_f__Intercept"     "r_ices_rect_f[39G4,Intercept]"
#> [21] "r_ices_rect_f[40G4,Intercept]" "r_ices_rect_f[40G5,Intercept]"
#> [23] "r_ices_rect_f[40G6,Intercept]" "r_ices_rect_f[40G7,Intercept]"
#> [25] "r_ices_rect_f[40G8,Intercept]" "r_ices_rect_f[41G4,Intercept]"
#> [27] "r_ices_rect_f[41G6,Intercept]" "r_ices_rect_f[41G7,Intercept]"
#> [29] "r_ices_rect_f[41G8,Intercept]" "r_ices_rect_f[41G9,Intercept]"
#> [31] "r_ices_rect_f[42G6,Intercept]" "r_ices_rect_f[42G7,Intercept]"
#> [33] "r_ices_rect_f[43G6,Intercept]" "r_ices_rect_f[43G7,Intercept]"
#> [35] "r_ices_rect_f[43G8,Intercept]" "r_ices_rect_f[43G9,Intercept]"
#> [37] "r_ices_rect_f[44G9,Intercept]" "lp__"                         
#> [39] "z_1[1,1]"                      "z_1[1,2]"                     
#> [41] "z_1[1,3]"                      "z_1[1,4]"                     
#> [43] "z_1[1,5]"                      "z_1[1,6]"                     
#> [45] "z_1[1,7]"                      "z_1[1,8]"                     
#> [47] "z_1[1,9]"                      "z_1[1,10]"                    
#> [49] "z_1[1,11]"                     "z_1[1,12]"                    
#> [51] "z_1[1,13]"                     "z_1[1,14]"                    
#> [53] "z_1[1,15]"                     "z_1[1,16]"                    
#> [55] "z_1[1,17]"                     "z_1[1,18]"
pal_diag <- rev(brewer.pal(n = 4, name = "Dark2"))

mcmc_trace(posterior,
           pars = c("b_year_f2015", "b_year_f2016", "b_year_f2017", "b_year_f2018",
                    "b_year_f2018", "b_year_f2019", "b_year_f2020",
                    "b_predfle_density_tot_sc", "b_predcod_density_tot_sc",
                    "b_quarter_f4", "b_phi_quarter_f4", 
                    "b_phi_year_f2015", "b_phi_year_f2016", "b_phi_year_f2017",
                    "b_phi_year_f2018", "b_phi_year_f2019", "b_phi_year_f2020", 
                    "b_phi_predfle_density_tot_sc", "b_phi_predcod_density_tot_sc", 
                    "sd_ices_rect_f__Intercept"),
                 facet_args = list(ncol = 2, strip.position = "left")) +
  theme(text = element_text(size = 12), strip.text = element_text(size = 4),
        legend.position = "top") +
  scale_color_manual(values = alpha(pal_diag, alpha = 0.6))
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.

ggsave("figures/supp/schoener_mcmc_trace.png", width = 6.5, height = 6.5, dpi = 600)

# Working with the posterior
posterior_m_schoen_beta <- m_schoen_beta %>% 
  gather_draws(`b_.*`, regex = TRUE)
  
ggplot(posterior_m_schoen_beta, aes(x = .value, y = fct_rev(.variable))) +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.7) +
  stat_halfeye(alpha = 0.5, .width = c(0.8, 0.95), point_interval = "median_hdi") +
  guides(fill = "none", slab_alpha = "none") +
  labs(x = "Coefficient", y = "Variable") +
  NULL

# Marginal effects of densities
m_schoen_beta_pred_cod <- m_schoen_beta %>% 
  epred_draws(newdata = tibble(predcod_density_tot_sc = 
                                 seq(min(schoener$predcod_density_tot_sc),
                                     max(schoener$predcod_density_tot_sc),
                                     length.out = 1000),
                               predfle_density_tot_sc = 0,
                               year_f = factor(2018),
                               quarter_f = factor(1)),
              re_formula = NA)

m_schoen_beta_pred_fle <- m_schoen_beta %>% 
  epred_draws(newdata = tibble(predfle_density_tot_sc = 
                                 seq(min(schoener$predfle_density_tot_sc),
                                     max(schoener$predfle_density_tot_sc),
                                     length.out = 1000),
                               predcod_density_tot_sc = 0,
                               year_f = factor(2018),
                               quarter_f = factor(1)),
              re_formula = NA)

p1 <- ggplot(m_schoen_beta_pred_fle, aes(x = predfle_density_tot_sc, y = .epred)) +
  stat_lineribbon() + 
  scale_fill_brewer(palette = "Purples") +
  theme_light(base_size = 14) +
  theme(legend.position = "bottom") +
  guides(fill = FALSE) +
  labs(x = "Scaled flounder density", y = NULL
       #, caption = "80% and 95% credible intervals shown in black"
       ) +
  NULL
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.

p2 <- ggplot(m_schoen_beta_pred_cod, aes(x = predcod_density_tot_sc, y = .epred)) +
  stat_lineribbon() + 
  scale_fill_brewer(palette = "Purples", name = "Credible interval") +
  theme_light(base_size = 14) +
  theme(legend.position = "bottom") +
  guides(fill = guide_legend(title.position = "top", title.hjust = 0.5)) +
  labs(x = "Scaled cod density", y = NULL
       , caption = "80% and 95% credible intervals shown in black"
       ) +
  NULL

p1/p2

ggsave("figures/schoener_dens_marginal.png", width = 6.5, height = 6.5, dpi = 600)

Beta model to size based schoener index

m_schoen_size_beta <- brm(
  bf(schoener2 ~ 0 + schoener_combination_f + year_f + quarter_f + (1|ices_rect_f),
     phi ~ 0 + schoener_combination_f + year_f + quarter_f),
  data = schoener_long2, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
  chains = 4, iter = 4000, cores = 4, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling
  
plot(m_schoen_size_beta)

loo_m_schoen_size_beta <- loo(m_schoen_size_beta, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.

Plot

# Evaluate fit and convergence etc.
# PP check
pp_check(m_schoen_size_beta, ndraws = 50) +
  theme_light(base_size = 20) + 
  scale_color_brewer(palette = "Dark2", name = "") +
  NULL
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.

ggsave("figures/supp/schoener_size_pp_check.png", width = 6.5, height = 6.5, dpi = 600)

# Chain convergence
posterior <- as.array(m_schoen_size_beta)
dimnames(posterior)
#> $iteration
#>    [1] "1"    "2"    "3"    "4"    "5"    "6"    "7"    "8"    "9"    "10"  
#>   [11] "11"   "12"   "13"   "14"   "15"   "16"   "17"   "18"   "19"   "20"  
#>   [21] "21"   "22"   "23"   "24"   "25"   "26"   "27"   "28"   "29"   "30"  
#>   [31] "31"   "32"   "33"   "34"   "35"   "36"   "37"   "38"   "39"   "40"  
#>   [41] "41"   "42"   "43"   "44"   "45"   "46"   "47"   "48"   "49"   "50"  
#>   [51] "51"   "52"   "53"   "54"   "55"   "56"   "57"   "58"   "59"   "60"  
#>   [61] "61"   "62"   "63"   "64"   "65"   "66"   "67"   "68"   "69"   "70"  
#>   [71] "71"   "72"   "73"   "74"   "75"   "76"   "77"   "78"   "79"   "80"  
#>   [81] "81"   "82"   "83"   "84"   "85"   "86"   "87"   "88"   "89"   "90"  
#>   [91] "91"   "92"   "93"   "94"   "95"   "96"   "97"   "98"   "99"   "100" 
#>  [101] "101"  "102"  "103"  "104"  "105"  "106"  "107"  "108"  "109"  "110" 
#>  [111] "111"  "112"  "113"  "114"  "115"  "116"  "117"  "118"  "119"  "120" 
#>  [121] "121"  "122"  "123"  "124"  "125"  "126"  "127"  "128"  "129"  "130" 
#>  [131] "131"  "132"  "133"  "134"  "135"  "136"  "137"  "138"  "139"  "140" 
#>  [141] "141"  "142"  "143"  "144"  "145"  "146"  "147"  "148"  "149"  "150" 
#>  [151] "151"  "152"  "153"  "154"  "155"  "156"  "157"  "158"  "159"  "160" 
#>  [161] "161"  "162"  "163"  "164"  "165"  "166"  "167"  "168"  "169"  "170" 
#>  [171] "171"  "172"  "173"  "174"  "175"  "176"  "177"  "178"  "179"  "180" 
#>  [181] "181"  "182"  "183"  "184"  "185"  "186"  "187"  "188"  "189"  "190" 
#>  [191] "191"  "192"  "193"  "194"  "195"  "196"  "197"  "198"  "199"  "200" 
#>  [201] "201"  "202"  "203"  "204"  "205"  "206"  "207"  "208"  "209"  "210" 
#>  [211] "211"  "212"  "213"  "214"  "215"  "216"  "217"  "218"  "219"  "220" 
#>  [221] "221"  "222"  "223"  "224"  "225"  "226"  "227"  "228"  "229"  "230" 
#>  [231] "231"  "232"  "233"  "234"  "235"  "236"  "237"  "238"  "239"  "240" 
#>  [241] "241"  "242"  "243"  "244"  "245"  "246"  "247"  "248"  "249"  "250" 
#>  [251] "251"  "252"  "253"  "254"  "255"  "256"  "257"  "258"  "259"  "260" 
#>  [261] "261"  "262"  "263"  "264"  "265"  "266"  "267"  "268"  "269"  "270" 
#>  [271] "271"  "272"  "273"  "274"  "275"  "276"  "277"  "278"  "279"  "280" 
#>  [281] "281"  "282"  "283"  "284"  "285"  "286"  "287"  "288"  "289"  "290" 
#>  [291] "291"  "292"  "293"  "294"  "295"  "296"  "297"  "298"  "299"  "300" 
#>  [301] "301"  "302"  "303"  "304"  "305"  "306"  "307"  "308"  "309"  "310" 
#>  [311] "311"  "312"  "313"  "314"  "315"  "316"  "317"  "318"  "319"  "320" 
#>  [321] "321"  "322"  "323"  "324"  "325"  "326"  "327"  "328"  "329"  "330" 
#>  [331] "331"  "332"  "333"  "334"  "335"  "336"  "337"  "338"  "339"  "340" 
#>  [341] "341"  "342"  "343"  "344"  "345"  "346"  "347"  "348"  "349"  "350" 
#>  [351] "351"  "352"  "353"  "354"  "355"  "356"  "357"  "358"  "359"  "360" 
#>  [361] "361"  "362"  "363"  "364"  "365"  "366"  "367"  "368"  "369"  "370" 
#>  [371] "371"  "372"  "373"  "374"  "375"  "376"  "377"  "378"  "379"  "380" 
#>  [381] "381"  "382"  "383"  "384"  "385"  "386"  "387"  "388"  "389"  "390" 
#>  [391] "391"  "392"  "393"  "394"  "395"  "396"  "397"  "398"  "399"  "400" 
#>  [401] "401"  "402"  "403"  "404"  "405"  "406"  "407"  "408"  "409"  "410" 
#>  [411] "411"  "412"  "413"  "414"  "415"  "416"  "417"  "418"  "419"  "420" 
#>  [421] "421"  "422"  "423"  "424"  "425"  "426"  "427"  "428"  "429"  "430" 
#>  [431] "431"  "432"  "433"  "434"  "435"  "436"  "437"  "438"  "439"  "440" 
#>  [441] "441"  "442"  "443"  "444"  "445"  "446"  "447"  "448"  "449"  "450" 
#>  [451] "451"  "452"  "453"  "454"  "455"  "456"  "457"  "458"  "459"  "460" 
#>  [461] "461"  "462"  "463"  "464"  "465"  "466"  "467"  "468"  "469"  "470" 
#>  [471] "471"  "472"  "473"  "474"  "475"  "476"  "477"  "478"  "479"  "480" 
#>  [481] "481"  "482"  "483"  "484"  "485"  "486"  "487"  "488"  "489"  "490" 
#>  [491] "491"  "492"  "493"  "494"  "495"  "496"  "497"  "498"  "499"  "500" 
#>  [501] "501"  "502"  "503"  "504"  "505"  "506"  "507"  "508"  "509"  "510" 
#>  [511] "511"  "512"  "513"  "514"  "515"  "516"  "517"  "518"  "519"  "520" 
#>  [521] "521"  "522"  "523"  "524"  "525"  "526"  "527"  "528"  "529"  "530" 
#>  [531] "531"  "532"  "533"  "534"  "535"  "536"  "537"  "538"  "539"  "540" 
#>  [541] "541"  "542"  "543"  "544"  "545"  "546"  "547"  "548"  "549"  "550" 
#>  [551] "551"  "552"  "553"  "554"  "555"  "556"  "557"  "558"  "559"  "560" 
#>  [561] "561"  "562"  "563"  "564"  "565"  "566"  "567"  "568"  "569"  "570" 
#>  [571] "571"  "572"  "573"  "574"  "575"  "576"  "577"  "578"  "579"  "580" 
#>  [581] "581"  "582"  "583"  "584"  "585"  "586"  "587"  "588"  "589"  "590" 
#>  [591] "591"  "592"  "593"  "594"  "595"  "596"  "597"  "598"  "599"  "600" 
#>  [601] "601"  "602"  "603"  "604"  "605"  "606"  "607"  "608"  "609"  "610" 
#>  [611] "611"  "612"  "613"  "614"  "615"  "616"  "617"  "618"  "619"  "620" 
#>  [621] "621"  "622"  "623"  "624"  "625"  "626"  "627"  "628"  "629"  "630" 
#>  [631] "631"  "632"  "633"  "634"  "635"  "636"  "637"  "638"  "639"  "640" 
#>  [641] "641"  "642"  "643"  "644"  "645"  "646"  "647"  "648"  "649"  "650" 
#>  [651] "651"  "652"  "653"  "654"  "655"  "656"  "657"  "658"  "659"  "660" 
#>  [661] "661"  "662"  "663"  "664"  "665"  "666"  "667"  "668"  "669"  "670" 
#>  [671] "671"  "672"  "673"  "674"  "675"  "676"  "677"  "678"  "679"  "680" 
#>  [681] "681"  "682"  "683"  "684"  "685"  "686"  "687"  "688"  "689"  "690" 
#>  [691] "691"  "692"  "693"  "694"  "695"  "696"  "697"  "698"  "699"  "700" 
#>  [701] "701"  "702"  "703"  "704"  "705"  "706"  "707"  "708"  "709"  "710" 
#>  [711] "711"  "712"  "713"  "714"  "715"  "716"  "717"  "718"  "719"  "720" 
#>  [721] "721"  "722"  "723"  "724"  "725"  "726"  "727"  "728"  "729"  "730" 
#>  [731] "731"  "732"  "733"  "734"  "735"  "736"  "737"  "738"  "739"  "740" 
#>  [741] "741"  "742"  "743"  "744"  "745"  "746"  "747"  "748"  "749"  "750" 
#>  [751] "751"  "752"  "753"  "754"  "755"  "756"  "757"  "758"  "759"  "760" 
#>  [761] "761"  "762"  "763"  "764"  "765"  "766"  "767"  "768"  "769"  "770" 
#>  [771] "771"  "772"  "773"  "774"  "775"  "776"  "777"  "778"  "779"  "780" 
#>  [781] "781"  "782"  "783"  "784"  "785"  "786"  "787"  "788"  "789"  "790" 
#>  [791] "791"  "792"  "793"  "794"  "795"  "796"  "797"  "798"  "799"  "800" 
#>  [801] "801"  "802"  "803"  "804"  "805"  "806"  "807"  "808"  "809"  "810" 
#>  [811] "811"  "812"  "813"  "814"  "815"  "816"  "817"  "818"  "819"  "820" 
#>  [821] "821"  "822"  "823"  "824"  "825"  "826"  "827"  "828"  "829"  "830" 
#>  [831] "831"  "832"  "833"  "834"  "835"  "836"  "837"  "838"  "839"  "840" 
#>  [841] "841"  "842"  "843"  "844"  "845"  "846"  "847"  "848"  "849"  "850" 
#>  [851] "851"  "852"  "853"  "854"  "855"  "856"  "857"  "858"  "859"  "860" 
#>  [861] "861"  "862"  "863"  "864"  "865"  "866"  "867"  "868"  "869"  "870" 
#>  [871] "871"  "872"  "873"  "874"  "875"  "876"  "877"  "878"  "879"  "880" 
#>  [881] "881"  "882"  "883"  "884"  "885"  "886"  "887"  "888"  "889"  "890" 
#>  [891] "891"  "892"  "893"  "894"  "895"  "896"  "897"  "898"  "899"  "900" 
#>  [901] "901"  "902"  "903"  "904"  "905"  "906"  "907"  "908"  "909"  "910" 
#>  [911] "911"  "912"  "913"  "914"  "915"  "916"  "917"  "918"  "919"  "920" 
#>  [921] "921"  "922"  "923"  "924"  "925"  "926"  "927"  "928"  "929"  "930" 
#>  [931] "931"  "932"  "933"  "934"  "935"  "936"  "937"  "938"  "939"  "940" 
#>  [941] "941"  "942"  "943"  "944"  "945"  "946"  "947"  "948"  "949"  "950" 
#>  [951] "951"  "952"  "953"  "954"  "955"  "956"  "957"  "958"  "959"  "960" 
#>  [961] "961"  "962"  "963"  "964"  "965"  "966"  "967"  "968"  "969"  "970" 
#>  [971] "971"  "972"  "973"  "974"  "975"  "976"  "977"  "978"  "979"  "980" 
#>  [981] "981"  "982"  "983"  "984"  "985"  "986"  "987"  "988"  "989"  "990" 
#>  [991] "991"  "992"  "993"  "994"  "995"  "996"  "997"  "998"  "999"  "1000"
#> [1001] "1001" "1002" "1003" "1004" "1005" "1006" "1007" "1008" "1009" "1010"
#> [1011] "1011" "1012" "1013" "1014" "1015" "1016" "1017" "1018" "1019" "1020"
#> [1021] "1021" "1022" "1023" "1024" "1025" "1026" "1027" "1028" "1029" "1030"
#> [1031] "1031" "1032" "1033" "1034" "1035" "1036" "1037" "1038" "1039" "1040"
#> [1041] "1041" "1042" "1043" "1044" "1045" "1046" "1047" "1048" "1049" "1050"
#> [1051] "1051" "1052" "1053" "1054" "1055" "1056" "1057" "1058" "1059" "1060"
#> [1061] "1061" "1062" "1063" "1064" "1065" "1066" "1067" "1068" "1069" "1070"
#> [1071] "1071" "1072" "1073" "1074" "1075" "1076" "1077" "1078" "1079" "1080"
#> [1081] "1081" "1082" "1083" "1084" "1085" "1086" "1087" "1088" "1089" "1090"
#> [1091] "1091" "1092" "1093" "1094" "1095" "1096" "1097" "1098" "1099" "1100"
#> [1101] "1101" "1102" "1103" "1104" "1105" "1106" "1107" "1108" "1109" "1110"
#> [1111] "1111" "1112" "1113" "1114" "1115" "1116" "1117" "1118" "1119" "1120"
#> [1121] "1121" "1122" "1123" "1124" "1125" "1126" "1127" "1128" "1129" "1130"
#> [1131] "1131" "1132" "1133" "1134" "1135" "1136" "1137" "1138" "1139" "1140"
#> [1141] "1141" "1142" "1143" "1144" "1145" "1146" "1147" "1148" "1149" "1150"
#> [1151] "1151" "1152" "1153" "1154" "1155" "1156" "1157" "1158" "1159" "1160"
#> [1161] "1161" "1162" "1163" "1164" "1165" "1166" "1167" "1168" "1169" "1170"
#> [1171] "1171" "1172" "1173" "1174" "1175" "1176" "1177" "1178" "1179" "1180"
#> [1181] "1181" "1182" "1183" "1184" "1185" "1186" "1187" "1188" "1189" "1190"
#> [1191] "1191" "1192" "1193" "1194" "1195" "1196" "1197" "1198" "1199" "1200"
#> [1201] "1201" "1202" "1203" "1204" "1205" "1206" "1207" "1208" "1209" "1210"
#> [1211] "1211" "1212" "1213" "1214" "1215" "1216" "1217" "1218" "1219" "1220"
#> [1221] "1221" "1222" "1223" "1224" "1225" "1226" "1227" "1228" "1229" "1230"
#> [1231] "1231" "1232" "1233" "1234" "1235" "1236" "1237" "1238" "1239" "1240"
#> [1241] "1241" "1242" "1243" "1244" "1245" "1246" "1247" "1248" "1249" "1250"
#> [1251] "1251" "1252" "1253" "1254" "1255" "1256" "1257" "1258" "1259" "1260"
#> [1261] "1261" "1262" "1263" "1264" "1265" "1266" "1267" "1268" "1269" "1270"
#> [1271] "1271" "1272" "1273" "1274" "1275" "1276" "1277" "1278" "1279" "1280"
#> [1281] "1281" "1282" "1283" "1284" "1285" "1286" "1287" "1288" "1289" "1290"
#> [1291] "1291" "1292" "1293" "1294" "1295" "1296" "1297" "1298" "1299" "1300"
#> [1301] "1301" "1302" "1303" "1304" "1305" "1306" "1307" "1308" "1309" "1310"
#> [1311] "1311" "1312" "1313" "1314" "1315" "1316" "1317" "1318" "1319" "1320"
#> [1321] "1321" "1322" "1323" "1324" "1325" "1326" "1327" "1328" "1329" "1330"
#> [1331] "1331" "1332" "1333" "1334" "1335" "1336" "1337" "1338" "1339" "1340"
#> [1341] "1341" "1342" "1343" "1344" "1345" "1346" "1347" "1348" "1349" "1350"
#> [1351] "1351" "1352" "1353" "1354" "1355" "1356" "1357" "1358" "1359" "1360"
#> [1361] "1361" "1362" "1363" "1364" "1365" "1366" "1367" "1368" "1369" "1370"
#> [1371] "1371" "1372" "1373" "1374" "1375" "1376" "1377" "1378" "1379" "1380"
#> [1381] "1381" "1382" "1383" "1384" "1385" "1386" "1387" "1388" "1389" "1390"
#> [1391] "1391" "1392" "1393" "1394" "1395" "1396" "1397" "1398" "1399" "1400"
#> [1401] "1401" "1402" "1403" "1404" "1405" "1406" "1407" "1408" "1409" "1410"
#> [1411] "1411" "1412" "1413" "1414" "1415" "1416" "1417" "1418" "1419" "1420"
#> [1421] "1421" "1422" "1423" "1424" "1425" "1426" "1427" "1428" "1429" "1430"
#> [1431] "1431" "1432" "1433" "1434" "1435" "1436" "1437" "1438" "1439" "1440"
#> [1441] "1441" "1442" "1443" "1444" "1445" "1446" "1447" "1448" "1449" "1450"
#> [1451] "1451" "1452" "1453" "1454" "1455" "1456" "1457" "1458" "1459" "1460"
#> [1461] "1461" "1462" "1463" "1464" "1465" "1466" "1467" "1468" "1469" "1470"
#> [1471] "1471" "1472" "1473" "1474" "1475" "1476" "1477" "1478" "1479" "1480"
#> [1481] "1481" "1482" "1483" "1484" "1485" "1486" "1487" "1488" "1489" "1490"
#> [1491] "1491" "1492" "1493" "1494" "1495" "1496" "1497" "1498" "1499" "1500"
#> [1501] "1501" "1502" "1503" "1504" "1505" "1506" "1507" "1508" "1509" "1510"
#> [1511] "1511" "1512" "1513" "1514" "1515" "1516" "1517" "1518" "1519" "1520"
#> [1521] "1521" "1522" "1523" "1524" "1525" "1526" "1527" "1528" "1529" "1530"
#> [1531] "1531" "1532" "1533" "1534" "1535" "1536" "1537" "1538" "1539" "1540"
#> [1541] "1541" "1542" "1543" "1544" "1545" "1546" "1547" "1548" "1549" "1550"
#> [1551] "1551" "1552" "1553" "1554" "1555" "1556" "1557" "1558" "1559" "1560"
#> [1561] "1561" "1562" "1563" "1564" "1565" "1566" "1567" "1568" "1569" "1570"
#> [1571] "1571" "1572" "1573" "1574" "1575" "1576" "1577" "1578" "1579" "1580"
#> [1581] "1581" "1582" "1583" "1584" "1585" "1586" "1587" "1588" "1589" "1590"
#> [1591] "1591" "1592" "1593" "1594" "1595" "1596" "1597" "1598" "1599" "1600"
#> [1601] "1601" "1602" "1603" "1604" "1605" "1606" "1607" "1608" "1609" "1610"
#> [1611] "1611" "1612" "1613" "1614" "1615" "1616" "1617" "1618" "1619" "1620"
#> [1621] "1621" "1622" "1623" "1624" "1625" "1626" "1627" "1628" "1629" "1630"
#> [1631] "1631" "1632" "1633" "1634" "1635" "1636" "1637" "1638" "1639" "1640"
#> [1641] "1641" "1642" "1643" "1644" "1645" "1646" "1647" "1648" "1649" "1650"
#> [1651] "1651" "1652" "1653" "1654" "1655" "1656" "1657" "1658" "1659" "1660"
#> [1661] "1661" "1662" "1663" "1664" "1665" "1666" "1667" "1668" "1669" "1670"
#> [1671] "1671" "1672" "1673" "1674" "1675" "1676" "1677" "1678" "1679" "1680"
#> [1681] "1681" "1682" "1683" "1684" "1685" "1686" "1687" "1688" "1689" "1690"
#> [1691] "1691" "1692" "1693" "1694" "1695" "1696" "1697" "1698" "1699" "1700"
#> [1701] "1701" "1702" "1703" "1704" "1705" "1706" "1707" "1708" "1709" "1710"
#> [1711] "1711" "1712" "1713" "1714" "1715" "1716" "1717" "1718" "1719" "1720"
#> [1721] "1721" "1722" "1723" "1724" "1725" "1726" "1727" "1728" "1729" "1730"
#> [1731] "1731" "1732" "1733" "1734" "1735" "1736" "1737" "1738" "1739" "1740"
#> [1741] "1741" "1742" "1743" "1744" "1745" "1746" "1747" "1748" "1749" "1750"
#> [1751] "1751" "1752" "1753" "1754" "1755" "1756" "1757" "1758" "1759" "1760"
#> [1761] "1761" "1762" "1763" "1764" "1765" "1766" "1767" "1768" "1769" "1770"
#> [1771] "1771" "1772" "1773" "1774" "1775" "1776" "1777" "1778" "1779" "1780"
#> [1781] "1781" "1782" "1783" "1784" "1785" "1786" "1787" "1788" "1789" "1790"
#> [1791] "1791" "1792" "1793" "1794" "1795" "1796" "1797" "1798" "1799" "1800"
#> [1801] "1801" "1802" "1803" "1804" "1805" "1806" "1807" "1808" "1809" "1810"
#> [1811] "1811" "1812" "1813" "1814" "1815" "1816" "1817" "1818" "1819" "1820"
#> [1821] "1821" "1822" "1823" "1824" "1825" "1826" "1827" "1828" "1829" "1830"
#> [1831] "1831" "1832" "1833" "1834" "1835" "1836" "1837" "1838" "1839" "1840"
#> [1841] "1841" "1842" "1843" "1844" "1845" "1846" "1847" "1848" "1849" "1850"
#> [1851] "1851" "1852" "1853" "1854" "1855" "1856" "1857" "1858" "1859" "1860"
#> [1861] "1861" "1862" "1863" "1864" "1865" "1866" "1867" "1868" "1869" "1870"
#> [1871] "1871" "1872" "1873" "1874" "1875" "1876" "1877" "1878" "1879" "1880"
#> [1881] "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889" "1890"
#> [1891] "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899" "1900"
#> [1901] "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909" "1910"
#> [1911] "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919" "1920"
#> [1921] "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929" "1930"
#> [1931] "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939" "1940"
#> [1941] "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949" "1950"
#> [1951] "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959" "1960"
#> [1961] "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969" "1970"
#> [1971] "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979" "1980"
#> [1981] "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989" "1990"
#> [1991] "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999" "2000"
#> 
#> $chain
#> [1] "1" "2" "3" "4"
#> 
#> $variable
#>  [1] "b_schoener_combination_fschoener_f_lc"     
#>  [2] "b_schoener_combination_fschoener_f_sc"     
#>  [3] "b_schoener_combination_fschoener_sc_lc"    
#>  [4] "b_year_f2016"                              
#>  [5] "b_year_f2017"                              
#>  [6] "b_year_f2018"                              
#>  [7] "b_year_f2019"                              
#>  [8] "b_year_f2020"                              
#>  [9] "b_quarter_f4"                              
#> [10] "b_phi_schoener_combination_fschoener_f_lc" 
#> [11] "b_phi_schoener_combination_fschoener_f_sc" 
#> [12] "b_phi_schoener_combination_fschoener_sc_lc"
#> [13] "b_phi_year_f2016"                          
#> [14] "b_phi_year_f2017"                          
#> [15] "b_phi_year_f2018"                          
#> [16] "b_phi_year_f2019"                          
#> [17] "b_phi_year_f2020"                          
#> [18] "b_phi_quarter_f4"                          
#> [19] "sd_ices_rect_f__Intercept"                 
#> [20] "r_ices_rect_f[39G4,Intercept]"             
#> [21] "r_ices_rect_f[40G4,Intercept]"             
#> [22] "r_ices_rect_f[40G5,Intercept]"             
#> [23] "r_ices_rect_f[40G6,Intercept]"             
#> [24] "r_ices_rect_f[40G7,Intercept]"             
#> [25] "r_ices_rect_f[40G8,Intercept]"             
#> [26] "r_ices_rect_f[41G4,Intercept]"             
#> [27] "r_ices_rect_f[41G6,Intercept]"             
#> [28] "r_ices_rect_f[41G7,Intercept]"             
#> [29] "r_ices_rect_f[41G8,Intercept]"             
#> [30] "r_ices_rect_f[41G9,Intercept]"             
#> [31] "r_ices_rect_f[42G6,Intercept]"             
#> [32] "r_ices_rect_f[42G7,Intercept]"             
#> [33] "r_ices_rect_f[43G6,Intercept]"             
#> [34] "r_ices_rect_f[43G7,Intercept]"             
#> [35] "r_ices_rect_f[43G8,Intercept]"             
#> [36] "r_ices_rect_f[43G9,Intercept]"             
#> [37] "r_ices_rect_f[44G9,Intercept]"             
#> [38] "lp__"                                      
#> [39] "z_1[1,1]"                                  
#> [40] "z_1[1,2]"                                  
#> [41] "z_1[1,3]"                                  
#> [42] "z_1[1,4]"                                  
#> [43] "z_1[1,5]"                                  
#> [44] "z_1[1,6]"                                  
#> [45] "z_1[1,7]"                                  
#> [46] "z_1[1,8]"                                  
#> [47] "z_1[1,9]"                                  
#> [48] "z_1[1,10]"                                 
#> [49] "z_1[1,11]"                                 
#> [50] "z_1[1,12]"                                 
#> [51] "z_1[1,13]"                                 
#> [52] "z_1[1,14]"                                 
#> [53] "z_1[1,15]"                                 
#> [54] "z_1[1,16]"                                 
#> [55] "z_1[1,17]"                                 
#> [56] "z_1[1,18]"
pal_diag <- rev(brewer.pal(n = 4, name = "Dark2"))

mcmc_trace(posterior,
           pars = c("b_schoener_combination_fschoener_f_lc",
                    "b_schoener_combination_fschoener_f_sc",
                    "b_schoener_combination_fschoener_sc_lc",
                    "b_year_f2016", "b_year_f2017", "b_year_f2018", "b_year_f2019",
                    "b_year_f2020", "b_quarter_f4", "b_phi_quarter_f4",
                    "b_phi_year_f2016", "b_phi_year_f2017", "b_phi_year_f2018",
                    "b_phi_year_f2019", "b_phi_year_f2020"),
                 facet_args = list(ncol = 2, strip.position = "left")) +
  theme(text = element_text(size = 12), strip.text = element_text(size = 4),
        legend.position = "top") +
  scale_color_manual(values = alpha(pal_diag, alpha = 0.6))
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.

ggsave("figures/supp/schoener_size_mcmc_trace.png", width = 6.5, height = 6.5, dpi = 600)

# Working with the posterior
posterior_m_schoen_size_beta <- m_schoen_size_beta %>% 
  gather_draws(`b_.*`, regex = TRUE)

ggplot(posterior_m_schoen_size_beta, aes(x = .value, y = fct_rev(.variable))) +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.7) +
  stat_halfeye(alpha = 0.5, .width = c(0.8, 0.95), point_interval = "median_hdi") +
  guides(fill = "none", slab_alpha = "none") +
  labs(x = "Coefficient", y = "Variable") +
  NULL

# Marginal effects of Schoener combination variable
m_schoen_size_beta_pred <- m_schoen_size_beta %>% 
  epred_draws(newdata = tibble(schoener_combination_f = c("schoener_f_sc", "schoener_f_lc", "schoener_sc_lc"),
                               year_f = factor(2018),
                               quarter_f = factor(1)),
              re_formula = NA) %>% 
  mutate(schoener_combination_f = ifelse(schoener_combination_f == "schoener_f_sc", 
                                         "Flounder-Cod (S)", schoener_combination_f),
         schoener_combination_f = ifelse(schoener_combination_f == "schoener_f_lc", 
                                         "Flounder-Cod (L)", schoener_combination_f),
         schoener_combination_f = ifelse(schoener_combination_f == "schoener_sc_lc", 
                                         "Cod (S)-Cod (L)", schoener_combination_f))
#> mutate (grouped): changed 24,000 values (100%) of 'schoener_combination_f' (0 new NA)

ggplot(m_schoen_size_beta_pred, aes(x = .epred, fill = schoener_combination_f)) +
  stat_halfeye(.width = c(0.75, 0.95), point_interval = "median_hdi", alpha = 0.5,
               position = position_dodge(width = 0.01)) +
  coord_cartesian(ylim = c(0, 0.3)) +
  scale_fill_brewer(palette = "Dark2", name = "Species-size combination") +
  theme_light(base_size = 14) +
  theme(legend.position = "top") +
  geom_vline(xintercept = 0.6, linetype = 2, color = "gray20") +
  guides(fill = guide_legend(title.position = "top", title.hjust = 0.5)) +
  labs(x = "Predicted Schoener's overlap index", y = NULL
       , caption = "80% and 95% credible intervals shown in black"
       ) +
  NULL

ggsave("figures/schoener_size_marginal.png", width = 6.5, height = 6.5, dpi = 600)

Beta model to levin index

# Cod model
m_levin_beta_cod <- brm(
  bf(levin_cod2 ~ 1,
     phi ~ 1),
  data = levin, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
  chains = 4, iter = 4000, cores = 4, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling

plot(m_levin_beta_cod)

loo_m_levin_beta_cod <- loo(m_levin_beta_cod)

# Flounder models
m_levin_beta_fle <- brm(
  bf(levin_fle2 ~ 1,
     phi ~ 1),
  data = levin, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
  chains = 4, iter = 4000, cores = 4, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> recompiling to avoid crashing R session
#> Start sampling

plot(m_levin_beta_fle)

loo_m_levin_beta_fle <- loo(m_levin_beta_fle)

# Pooled model of Levin index and species as covariate
colnames(levin)
#>  [1] "year"       "quarter"    "area"       "levin_cod"  "levin_fle" 
#>  [6] "levin_cod2" "levin_fle2" "year_f"     "quarter_f"  "area_f"
levin2 <- levin %>% pivot_longer(4:5, names_to = "species", values_to = "levins")
#> pivot_longer: reorganized (levin_cod, levin_fle) into (species, levins) [was 22x10, now 44x10]
colnames(levin2)
#>  [1] "year"       "quarter"    "area"       "levin_cod2" "levin_fle2"
#>  [6] "year_f"     "quarter_f"  "area_f"     "species"    "levins"

m_levin_beta <- brm(
  bf(levins ~ 0 + species,
     phi ~ 0 + species),
  data = levin2, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
  chains = 4, iter = 4000, cores = 4, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling

plot(m_levin_beta)

loo_m_levin_beta <- loo(m_levin_beta)

Plot

# Evaluate fit and convergence etc.
# PP check
pp_check(m_levin_beta, ndraws = 50) +
  theme_light(base_size = 20) + 
  scale_color_brewer(palette = "Dark2", name = "") +
  NULL
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.

ggsave("figures/supp/levin_pp_check.png", width = 6.5, height = 6.5, dpi = 600)

# Chain convergence
posterior <- as.array(m_levin_beta)

dimnames(posterior)
#> $iteration
#>    [1] "1"    "2"    "3"    "4"    "5"    "6"    "7"    "8"    "9"    "10"  
#>   [11] "11"   "12"   "13"   "14"   "15"   "16"   "17"   "18"   "19"   "20"  
#>   [21] "21"   "22"   "23"   "24"   "25"   "26"   "27"   "28"   "29"   "30"  
#>   [31] "31"   "32"   "33"   "34"   "35"   "36"   "37"   "38"   "39"   "40"  
#>   [41] "41"   "42"   "43"   "44"   "45"   "46"   "47"   "48"   "49"   "50"  
#>   [51] "51"   "52"   "53"   "54"   "55"   "56"   "57"   "58"   "59"   "60"  
#>   [61] "61"   "62"   "63"   "64"   "65"   "66"   "67"   "68"   "69"   "70"  
#>   [71] "71"   "72"   "73"   "74"   "75"   "76"   "77"   "78"   "79"   "80"  
#>   [81] "81"   "82"   "83"   "84"   "85"   "86"   "87"   "88"   "89"   "90"  
#>   [91] "91"   "92"   "93"   "94"   "95"   "96"   "97"   "98"   "99"   "100" 
#>  [101] "101"  "102"  "103"  "104"  "105"  "106"  "107"  "108"  "109"  "110" 
#>  [111] "111"  "112"  "113"  "114"  "115"  "116"  "117"  "118"  "119"  "120" 
#>  [121] "121"  "122"  "123"  "124"  "125"  "126"  "127"  "128"  "129"  "130" 
#>  [131] "131"  "132"  "133"  "134"  "135"  "136"  "137"  "138"  "139"  "140" 
#>  [141] "141"  "142"  "143"  "144"  "145"  "146"  "147"  "148"  "149"  "150" 
#>  [151] "151"  "152"  "153"  "154"  "155"  "156"  "157"  "158"  "159"  "160" 
#>  [161] "161"  "162"  "163"  "164"  "165"  "166"  "167"  "168"  "169"  "170" 
#>  [171] "171"  "172"  "173"  "174"  "175"  "176"  "177"  "178"  "179"  "180" 
#>  [181] "181"  "182"  "183"  "184"  "185"  "186"  "187"  "188"  "189"  "190" 
#>  [191] "191"  "192"  "193"  "194"  "195"  "196"  "197"  "198"  "199"  "200" 
#>  [201] "201"  "202"  "203"  "204"  "205"  "206"  "207"  "208"  "209"  "210" 
#>  [211] "211"  "212"  "213"  "214"  "215"  "216"  "217"  "218"  "219"  "220" 
#>  [221] "221"  "222"  "223"  "224"  "225"  "226"  "227"  "228"  "229"  "230" 
#>  [231] "231"  "232"  "233"  "234"  "235"  "236"  "237"  "238"  "239"  "240" 
#>  [241] "241"  "242"  "243"  "244"  "245"  "246"  "247"  "248"  "249"  "250" 
#>  [251] "251"  "252"  "253"  "254"  "255"  "256"  "257"  "258"  "259"  "260" 
#>  [261] "261"  "262"  "263"  "264"  "265"  "266"  "267"  "268"  "269"  "270" 
#>  [271] "271"  "272"  "273"  "274"  "275"  "276"  "277"  "278"  "279"  "280" 
#>  [281] "281"  "282"  "283"  "284"  "285"  "286"  "287"  "288"  "289"  "290" 
#>  [291] "291"  "292"  "293"  "294"  "295"  "296"  "297"  "298"  "299"  "300" 
#>  [301] "301"  "302"  "303"  "304"  "305"  "306"  "307"  "308"  "309"  "310" 
#>  [311] "311"  "312"  "313"  "314"  "315"  "316"  "317"  "318"  "319"  "320" 
#>  [321] "321"  "322"  "323"  "324"  "325"  "326"  "327"  "328"  "329"  "330" 
#>  [331] "331"  "332"  "333"  "334"  "335"  "336"  "337"  "338"  "339"  "340" 
#>  [341] "341"  "342"  "343"  "344"  "345"  "346"  "347"  "348"  "349"  "350" 
#>  [351] "351"  "352"  "353"  "354"  "355"  "356"  "357"  "358"  "359"  "360" 
#>  [361] "361"  "362"  "363"  "364"  "365"  "366"  "367"  "368"  "369"  "370" 
#>  [371] "371"  "372"  "373"  "374"  "375"  "376"  "377"  "378"  "379"  "380" 
#>  [381] "381"  "382"  "383"  "384"  "385"  "386"  "387"  "388"  "389"  "390" 
#>  [391] "391"  "392"  "393"  "394"  "395"  "396"  "397"  "398"  "399"  "400" 
#>  [401] "401"  "402"  "403"  "404"  "405"  "406"  "407"  "408"  "409"  "410" 
#>  [411] "411"  "412"  "413"  "414"  "415"  "416"  "417"  "418"  "419"  "420" 
#>  [421] "421"  "422"  "423"  "424"  "425"  "426"  "427"  "428"  "429"  "430" 
#>  [431] "431"  "432"  "433"  "434"  "435"  "436"  "437"  "438"  "439"  "440" 
#>  [441] "441"  "442"  "443"  "444"  "445"  "446"  "447"  "448"  "449"  "450" 
#>  [451] "451"  "452"  "453"  "454"  "455"  "456"  "457"  "458"  "459"  "460" 
#>  [461] "461"  "462"  "463"  "464"  "465"  "466"  "467"  "468"  "469"  "470" 
#>  [471] "471"  "472"  "473"  "474"  "475"  "476"  "477"  "478"  "479"  "480" 
#>  [481] "481"  "482"  "483"  "484"  "485"  "486"  "487"  "488"  "489"  "490" 
#>  [491] "491"  "492"  "493"  "494"  "495"  "496"  "497"  "498"  "499"  "500" 
#>  [501] "501"  "502"  "503"  "504"  "505"  "506"  "507"  "508"  "509"  "510" 
#>  [511] "511"  "512"  "513"  "514"  "515"  "516"  "517"  "518"  "519"  "520" 
#>  [521] "521"  "522"  "523"  "524"  "525"  "526"  "527"  "528"  "529"  "530" 
#>  [531] "531"  "532"  "533"  "534"  "535"  "536"  "537"  "538"  "539"  "540" 
#>  [541] "541"  "542"  "543"  "544"  "545"  "546"  "547"  "548"  "549"  "550" 
#>  [551] "551"  "552"  "553"  "554"  "555"  "556"  "557"  "558"  "559"  "560" 
#>  [561] "561"  "562"  "563"  "564"  "565"  "566"  "567"  "568"  "569"  "570" 
#>  [571] "571"  "572"  "573"  "574"  "575"  "576"  "577"  "578"  "579"  "580" 
#>  [581] "581"  "582"  "583"  "584"  "585"  "586"  "587"  "588"  "589"  "590" 
#>  [591] "591"  "592"  "593"  "594"  "595"  "596"  "597"  "598"  "599"  "600" 
#>  [601] "601"  "602"  "603"  "604"  "605"  "606"  "607"  "608"  "609"  "610" 
#>  [611] "611"  "612"  "613"  "614"  "615"  "616"  "617"  "618"  "619"  "620" 
#>  [621] "621"  "622"  "623"  "624"  "625"  "626"  "627"  "628"  "629"  "630" 
#>  [631] "631"  "632"  "633"  "634"  "635"  "636"  "637"  "638"  "639"  "640" 
#>  [641] "641"  "642"  "643"  "644"  "645"  "646"  "647"  "648"  "649"  "650" 
#>  [651] "651"  "652"  "653"  "654"  "655"  "656"  "657"  "658"  "659"  "660" 
#>  [661] "661"  "662"  "663"  "664"  "665"  "666"  "667"  "668"  "669"  "670" 
#>  [671] "671"  "672"  "673"  "674"  "675"  "676"  "677"  "678"  "679"  "680" 
#>  [681] "681"  "682"  "683"  "684"  "685"  "686"  "687"  "688"  "689"  "690" 
#>  [691] "691"  "692"  "693"  "694"  "695"  "696"  "697"  "698"  "699"  "700" 
#>  [701] "701"  "702"  "703"  "704"  "705"  "706"  "707"  "708"  "709"  "710" 
#>  [711] "711"  "712"  "713"  "714"  "715"  "716"  "717"  "718"  "719"  "720" 
#>  [721] "721"  "722"  "723"  "724"  "725"  "726"  "727"  "728"  "729"  "730" 
#>  [731] "731"  "732"  "733"  "734"  "735"  "736"  "737"  "738"  "739"  "740" 
#>  [741] "741"  "742"  "743"  "744"  "745"  "746"  "747"  "748"  "749"  "750" 
#>  [751] "751"  "752"  "753"  "754"  "755"  "756"  "757"  "758"  "759"  "760" 
#>  [761] "761"  "762"  "763"  "764"  "765"  "766"  "767"  "768"  "769"  "770" 
#>  [771] "771"  "772"  "773"  "774"  "775"  "776"  "777"  "778"  "779"  "780" 
#>  [781] "781"  "782"  "783"  "784"  "785"  "786"  "787"  "788"  "789"  "790" 
#>  [791] "791"  "792"  "793"  "794"  "795"  "796"  "797"  "798"  "799"  "800" 
#>  [801] "801"  "802"  "803"  "804"  "805"  "806"  "807"  "808"  "809"  "810" 
#>  [811] "811"  "812"  "813"  "814"  "815"  "816"  "817"  "818"  "819"  "820" 
#>  [821] "821"  "822"  "823"  "824"  "825"  "826"  "827"  "828"  "829"  "830" 
#>  [831] "831"  "832"  "833"  "834"  "835"  "836"  "837"  "838"  "839"  "840" 
#>  [841] "841"  "842"  "843"  "844"  "845"  "846"  "847"  "848"  "849"  "850" 
#>  [851] "851"  "852"  "853"  "854"  "855"  "856"  "857"  "858"  "859"  "860" 
#>  [861] "861"  "862"  "863"  "864"  "865"  "866"  "867"  "868"  "869"  "870" 
#>  [871] "871"  "872"  "873"  "874"  "875"  "876"  "877"  "878"  "879"  "880" 
#>  [881] "881"  "882"  "883"  "884"  "885"  "886"  "887"  "888"  "889"  "890" 
#>  [891] "891"  "892"  "893"  "894"  "895"  "896"  "897"  "898"  "899"  "900" 
#>  [901] "901"  "902"  "903"  "904"  "905"  "906"  "907"  "908"  "909"  "910" 
#>  [911] "911"  "912"  "913"  "914"  "915"  "916"  "917"  "918"  "919"  "920" 
#>  [921] "921"  "922"  "923"  "924"  "925"  "926"  "927"  "928"  "929"  "930" 
#>  [931] "931"  "932"  "933"  "934"  "935"  "936"  "937"  "938"  "939"  "940" 
#>  [941] "941"  "942"  "943"  "944"  "945"  "946"  "947"  "948"  "949"  "950" 
#>  [951] "951"  "952"  "953"  "954"  "955"  "956"  "957"  "958"  "959"  "960" 
#>  [961] "961"  "962"  "963"  "964"  "965"  "966"  "967"  "968"  "969"  "970" 
#>  [971] "971"  "972"  "973"  "974"  "975"  "976"  "977"  "978"  "979"  "980" 
#>  [981] "981"  "982"  "983"  "984"  "985"  "986"  "987"  "988"  "989"  "990" 
#>  [991] "991"  "992"  "993"  "994"  "995"  "996"  "997"  "998"  "999"  "1000"
#> [1001] "1001" "1002" "1003" "1004" "1005" "1006" "1007" "1008" "1009" "1010"
#> [1011] "1011" "1012" "1013" "1014" "1015" "1016" "1017" "1018" "1019" "1020"
#> [1021] "1021" "1022" "1023" "1024" "1025" "1026" "1027" "1028" "1029" "1030"
#> [1031] "1031" "1032" "1033" "1034" "1035" "1036" "1037" "1038" "1039" "1040"
#> [1041] "1041" "1042" "1043" "1044" "1045" "1046" "1047" "1048" "1049" "1050"
#> [1051] "1051" "1052" "1053" "1054" "1055" "1056" "1057" "1058" "1059" "1060"
#> [1061] "1061" "1062" "1063" "1064" "1065" "1066" "1067" "1068" "1069" "1070"
#> [1071] "1071" "1072" "1073" "1074" "1075" "1076" "1077" "1078" "1079" "1080"
#> [1081] "1081" "1082" "1083" "1084" "1085" "1086" "1087" "1088" "1089" "1090"
#> [1091] "1091" "1092" "1093" "1094" "1095" "1096" "1097" "1098" "1099" "1100"
#> [1101] "1101" "1102" "1103" "1104" "1105" "1106" "1107" "1108" "1109" "1110"
#> [1111] "1111" "1112" "1113" "1114" "1115" "1116" "1117" "1118" "1119" "1120"
#> [1121] "1121" "1122" "1123" "1124" "1125" "1126" "1127" "1128" "1129" "1130"
#> [1131] "1131" "1132" "1133" "1134" "1135" "1136" "1137" "1138" "1139" "1140"
#> [1141] "1141" "1142" "1143" "1144" "1145" "1146" "1147" "1148" "1149" "1150"
#> [1151] "1151" "1152" "1153" "1154" "1155" "1156" "1157" "1158" "1159" "1160"
#> [1161] "1161" "1162" "1163" "1164" "1165" "1166" "1167" "1168" "1169" "1170"
#> [1171] "1171" "1172" "1173" "1174" "1175" "1176" "1177" "1178" "1179" "1180"
#> [1181] "1181" "1182" "1183" "1184" "1185" "1186" "1187" "1188" "1189" "1190"
#> [1191] "1191" "1192" "1193" "1194" "1195" "1196" "1197" "1198" "1199" "1200"
#> [1201] "1201" "1202" "1203" "1204" "1205" "1206" "1207" "1208" "1209" "1210"
#> [1211] "1211" "1212" "1213" "1214" "1215" "1216" "1217" "1218" "1219" "1220"
#> [1221] "1221" "1222" "1223" "1224" "1225" "1226" "1227" "1228" "1229" "1230"
#> [1231] "1231" "1232" "1233" "1234" "1235" "1236" "1237" "1238" "1239" "1240"
#> [1241] "1241" "1242" "1243" "1244" "1245" "1246" "1247" "1248" "1249" "1250"
#> [1251] "1251" "1252" "1253" "1254" "1255" "1256" "1257" "1258" "1259" "1260"
#> [1261] "1261" "1262" "1263" "1264" "1265" "1266" "1267" "1268" "1269" "1270"
#> [1271] "1271" "1272" "1273" "1274" "1275" "1276" "1277" "1278" "1279" "1280"
#> [1281] "1281" "1282" "1283" "1284" "1285" "1286" "1287" "1288" "1289" "1290"
#> [1291] "1291" "1292" "1293" "1294" "1295" "1296" "1297" "1298" "1299" "1300"
#> [1301] "1301" "1302" "1303" "1304" "1305" "1306" "1307" "1308" "1309" "1310"
#> [1311] "1311" "1312" "1313" "1314" "1315" "1316" "1317" "1318" "1319" "1320"
#> [1321] "1321" "1322" "1323" "1324" "1325" "1326" "1327" "1328" "1329" "1330"
#> [1331] "1331" "1332" "1333" "1334" "1335" "1336" "1337" "1338" "1339" "1340"
#> [1341] "1341" "1342" "1343" "1344" "1345" "1346" "1347" "1348" "1349" "1350"
#> [1351] "1351" "1352" "1353" "1354" "1355" "1356" "1357" "1358" "1359" "1360"
#> [1361] "1361" "1362" "1363" "1364" "1365" "1366" "1367" "1368" "1369" "1370"
#> [1371] "1371" "1372" "1373" "1374" "1375" "1376" "1377" "1378" "1379" "1380"
#> [1381] "1381" "1382" "1383" "1384" "1385" "1386" "1387" "1388" "1389" "1390"
#> [1391] "1391" "1392" "1393" "1394" "1395" "1396" "1397" "1398" "1399" "1400"
#> [1401] "1401" "1402" "1403" "1404" "1405" "1406" "1407" "1408" "1409" "1410"
#> [1411] "1411" "1412" "1413" "1414" "1415" "1416" "1417" "1418" "1419" "1420"
#> [1421] "1421" "1422" "1423" "1424" "1425" "1426" "1427" "1428" "1429" "1430"
#> [1431] "1431" "1432" "1433" "1434" "1435" "1436" "1437" "1438" "1439" "1440"
#> [1441] "1441" "1442" "1443" "1444" "1445" "1446" "1447" "1448" "1449" "1450"
#> [1451] "1451" "1452" "1453" "1454" "1455" "1456" "1457" "1458" "1459" "1460"
#> [1461] "1461" "1462" "1463" "1464" "1465" "1466" "1467" "1468" "1469" "1470"
#> [1471] "1471" "1472" "1473" "1474" "1475" "1476" "1477" "1478" "1479" "1480"
#> [1481] "1481" "1482" "1483" "1484" "1485" "1486" "1487" "1488" "1489" "1490"
#> [1491] "1491" "1492" "1493" "1494" "1495" "1496" "1497" "1498" "1499" "1500"
#> [1501] "1501" "1502" "1503" "1504" "1505" "1506" "1507" "1508" "1509" "1510"
#> [1511] "1511" "1512" "1513" "1514" "1515" "1516" "1517" "1518" "1519" "1520"
#> [1521] "1521" "1522" "1523" "1524" "1525" "1526" "1527" "1528" "1529" "1530"
#> [1531] "1531" "1532" "1533" "1534" "1535" "1536" "1537" "1538" "1539" "1540"
#> [1541] "1541" "1542" "1543" "1544" "1545" "1546" "1547" "1548" "1549" "1550"
#> [1551] "1551" "1552" "1553" "1554" "1555" "1556" "1557" "1558" "1559" "1560"
#> [1561] "1561" "1562" "1563" "1564" "1565" "1566" "1567" "1568" "1569" "1570"
#> [1571] "1571" "1572" "1573" "1574" "1575" "1576" "1577" "1578" "1579" "1580"
#> [1581] "1581" "1582" "1583" "1584" "1585" "1586" "1587" "1588" "1589" "1590"
#> [1591] "1591" "1592" "1593" "1594" "1595" "1596" "1597" "1598" "1599" "1600"
#> [1601] "1601" "1602" "1603" "1604" "1605" "1606" "1607" "1608" "1609" "1610"
#> [1611] "1611" "1612" "1613" "1614" "1615" "1616" "1617" "1618" "1619" "1620"
#> [1621] "1621" "1622" "1623" "1624" "1625" "1626" "1627" "1628" "1629" "1630"
#> [1631] "1631" "1632" "1633" "1634" "1635" "1636" "1637" "1638" "1639" "1640"
#> [1641] "1641" "1642" "1643" "1644" "1645" "1646" "1647" "1648" "1649" "1650"
#> [1651] "1651" "1652" "1653" "1654" "1655" "1656" "1657" "1658" "1659" "1660"
#> [1661] "1661" "1662" "1663" "1664" "1665" "1666" "1667" "1668" "1669" "1670"
#> [1671] "1671" "1672" "1673" "1674" "1675" "1676" "1677" "1678" "1679" "1680"
#> [1681] "1681" "1682" "1683" "1684" "1685" "1686" "1687" "1688" "1689" "1690"
#> [1691] "1691" "1692" "1693" "1694" "1695" "1696" "1697" "1698" "1699" "1700"
#> [1701] "1701" "1702" "1703" "1704" "1705" "1706" "1707" "1708" "1709" "1710"
#> [1711] "1711" "1712" "1713" "1714" "1715" "1716" "1717" "1718" "1719" "1720"
#> [1721] "1721" "1722" "1723" "1724" "1725" "1726" "1727" "1728" "1729" "1730"
#> [1731] "1731" "1732" "1733" "1734" "1735" "1736" "1737" "1738" "1739" "1740"
#> [1741] "1741" "1742" "1743" "1744" "1745" "1746" "1747" "1748" "1749" "1750"
#> [1751] "1751" "1752" "1753" "1754" "1755" "1756" "1757" "1758" "1759" "1760"
#> [1761] "1761" "1762" "1763" "1764" "1765" "1766" "1767" "1768" "1769" "1770"
#> [1771] "1771" "1772" "1773" "1774" "1775" "1776" "1777" "1778" "1779" "1780"
#> [1781] "1781" "1782" "1783" "1784" "1785" "1786" "1787" "1788" "1789" "1790"
#> [1791] "1791" "1792" "1793" "1794" "1795" "1796" "1797" "1798" "1799" "1800"
#> [1801] "1801" "1802" "1803" "1804" "1805" "1806" "1807" "1808" "1809" "1810"
#> [1811] "1811" "1812" "1813" "1814" "1815" "1816" "1817" "1818" "1819" "1820"
#> [1821] "1821" "1822" "1823" "1824" "1825" "1826" "1827" "1828" "1829" "1830"
#> [1831] "1831" "1832" "1833" "1834" "1835" "1836" "1837" "1838" "1839" "1840"
#> [1841] "1841" "1842" "1843" "1844" "1845" "1846" "1847" "1848" "1849" "1850"
#> [1851] "1851" "1852" "1853" "1854" "1855" "1856" "1857" "1858" "1859" "1860"
#> [1861] "1861" "1862" "1863" "1864" "1865" "1866" "1867" "1868" "1869" "1870"
#> [1871] "1871" "1872" "1873" "1874" "1875" "1876" "1877" "1878" "1879" "1880"
#> [1881] "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889" "1890"
#> [1891] "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899" "1900"
#> [1901] "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909" "1910"
#> [1911] "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919" "1920"
#> [1921] "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929" "1930"
#> [1931] "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939" "1940"
#> [1941] "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949" "1950"
#> [1951] "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959" "1960"
#> [1961] "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969" "1970"
#> [1971] "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979" "1980"
#> [1981] "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989" "1990"
#> [1991] "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999" "2000"
#> 
#> $chain
#> [1] "1" "2" "3" "4"
#> 
#> $variable
#> [1] "b_specieslevin_cod"     "b_specieslevin_fle"     "b_phi_specieslevin_cod"
#> [4] "b_phi_specieslevin_fle" "lp__"
pal_diag <- rev(brewer.pal(n = 4, name = "Dark2"))

mcmc_trace(posterior,
           pars = c("b_specieslevin_cod", "b_specieslevin_fle",
                    "b_phi_specieslevin_cod", "b_phi_specieslevin_fle"),
                 facet_args = list(ncol = 1, strip.position = "left")) +
  theme(text = element_text(size = 12), strip.text = element_text(size = 10),
        legend.position = "top") +
  scale_color_manual(values = alpha(pal_diag, alpha = 0.6))
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.

ggsave("figures/supp/levin_mcmc_trace.png", width = 6.5, height = 6.5, dpi = 600)

# Marginal effects of Levin's species covariate
beta_levin_pred <- m_levin_beta %>% 
    epred_draws(newdata = tibble(species = c("levin_cod", "levin_fle")),
                re_formula = NA) %>% 
  mutate(species = ifelse(species == "levin_fle", "Flounder", "Cod"))
#> mutate (grouped): changed 16,000 values (100%) of 'species' (0 new NA)

p1 <- ggplot(beta_levin_pred, aes(x = .epred, fill = species)) +
  stat_halfeye(.width = c(0.8, 0.95), point_interval = "median_hdi", alpha = 0.5,
               adjust = 2, position = position_dodge(width = 0.01)) +
  coord_cartesian(ylim = c(0, 0.45)) +
  scale_fill_brewer(palette = "Dark2", name = "Species") +
  theme_light(base_size = 16) +
  theme(legend.position = "top") +
  guides(fill = guide_legend(title.position = "top", title.hjust = 0.5)) +
  labs(x = "Predicted Levin's diversity index", y = NULL
       #, caption = "80% and 95% credible intervals shown in black"
       ) +
  NULL

#ggsave("figures/levin_marginal.png", width = 6.5, height = 6.5, dpi = 600)

# Now plot the difference!
beta_levin_pred_cod <- m_levin_beta %>% 
  epred_draws(newdata = tibble(species = "levin_cod"),
              re_formula = NA) %>% 
  rename(".epred_cod" = ".epred") %>% 
  ungroup() %>% 
  dplyr::select(-species)
#> rename: renamed one variable (.epred_cod)
#> ungroup: no grouping variables

beta_levin_pred_fle <- m_levin_beta %>% 
  epred_draws(newdata = tibble(species = "levin_fle"),
              re_formula = NA) %>% 
  rename(".epred_fle" = ".epred") %>% 
  ungroup() %>% 
  dplyr::select(-species)
#> rename: renamed one variable (.epred_fle)
#> ungroup: no grouping variables

beta_levin_pred_wide <- left_join(beta_levin_pred_cod, beta_levin_pred_fle) %>% 
  mutate(.epred_diff = .epred_cod - .epred_fle)
#> Joining, by = c(".row", ".chain", ".iteration", ".draw")
#> left_join: added one column (.epred_fle)
#>            > rows only in x       0
#>            > rows only in y  (    0)
#>            > matched rows     8,000
#>            >                 =======
#>            > rows total       8,000
#> mutate: new variable '.epred_diff' (double) with 8,000 unique values and 0% NA

p2 <- ggplot(beta_levin_pred_wide, aes(x = .epred_diff)) +
  stat_dotsinterval(quantiles = 100) + 
  geom_vline(xintercept = 0, linetype = 2, color = "tomato", size = 1) + 
  theme_light(base_size = 16) +
  guides(fill = guide_legend(title.position = "top", title.hjust = 0.5)) +
  labs(x = "Predicted difference in Levin's diversity index", y = NULL
       , caption = "80% and 95% credible intervals shown in black"
       ) +
  NULL

p1 / p2 #+ plot_annotation(tag_levels = "A")

ggsave("figures/levin_species_combined.png", width = 6.5, height = 6.5, dpi = 600)

TEST: Calculate Levins’s index per 5 cm length groups!

# This is just very exploratory...

fle_prey_long3 <- fle_prey %>%
  mutate(size_class = cut(pred_cm, breaks = c(seq(0, 100, 5)))) %>% 
  pivot_longer(15:30) %>% 
  group_by(name, year, quarter, area, size_class) %>%
  summarise(fle_stomach_content = sum(value)) %>% 
  arrange(name, size_class, year, area, quarter) %>%
  ungroup() %>% 
  mutate(id = paste(year, quarter, area, size_class, sep = "_")) %>% 
  group_by(id) %>%
  mutate(fle_stomach_content_tot = sum(fle_stomach_content),
         fle_stomach_content_prop = fle_stomach_content / fle_stomach_content_tot) %>% 
  ungroup()
#> mutate: new variable 'size_class' (factor) with 8 unique values and 0% NA
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 2180x33, now 34880x19]
#> group_by: 5 grouping variables (name, year, quarter, area, size_class)
#> summarise: now 1,952 rows and 6 columns, 4 group variables remaining (name, year, quarter, area)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 122 unique values and 0% NA
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'fle_stomach_content_tot' (double) with 121 unique values and 0% NA
#>                   new variable 'fle_stomach_content_prop' (double) with 676 unique values and 2% NA
#> ungroup: no grouping variables

# Now cod
cod_prey_long3 <- cod_prey %>%
  mutate(size_class = cut(pred_cm, breaks = c(seq(0, 100, 5)))) %>% 
  pivot_longer(15:30) %>% 
  group_by(name, year, quarter, area, size_class) %>%
  summarise(cod_stomach_content = sum(value)) %>% 
  arrange(name, size_class, year, area, quarter) %>%
  ungroup() %>% 
  mutate(id = paste(year, quarter, area, size_class, sep = "_")) %>% 
  group_by(id) %>%
  mutate(cod_stomach_content_tot = sum(cod_stomach_content),
         cod_stomach_content_prop = cod_stomach_content / cod_stomach_content_tot) %>% 
  ungroup()
#> mutate: new variable 'size_class' (factor) with 16 unique values and 0% NA
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x33, now 57488x19]
#> group_by: 5 grouping variables (name, year, quarter, area, size_class)
#> summarise: now 4,080 rows and 6 columns, 4 group variables remaining (name, year, quarter, area)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 255 unique values and 0% NA
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'cod_stomach_content_tot' (double) with 242 unique values and 0% NA
#>                   new variable 'cod_stomach_content_prop' (double) with 1,346 unique values and 1% NA
#> ungroup: no grouping variables

# Size based Levin's index
levin_cod <- cod_prey_long3 %>% 
  drop_na(name) %>% 
  drop_na(cod_stomach_content_prop) %>% 
  group_by(id, year, quarter, area, size_class) %>% 
  summarise(levin = ((1/(number_of_prey$n-1)) * (((1/sum(cod_stomach_content_prop^2))) - 1))) %>% 
  ungroup() %>% 
  mutate(size = as.integer(stringr::str_extract(size_class, "\\d+")))
#> drop_na: no rows removed
#> drop_na: removed 48 rows (1%), 4,032 rows remaining
#> group_by: 5 grouping variables (id, year, quarter, area, size_class)
#> summarise: now 252 rows and 6 columns, 4 group variables remaining (id, year, quarter, area)
#> ungroup: no grouping variables
#> mutate: new variable 'size' (integer) with 16 unique values and 0% NA

levin_fle <- fle_prey_long3 %>% 
  drop_na(name) %>% 
  drop_na(fle_stomach_content_prop) %>% 
  group_by(id, year, quarter, area, size_class) %>% 
  summarise(levin = ((1/(number_of_prey$n-1)) * (((1/sum(fle_stomach_content_prop^2))) - 1))) %>% 
  ungroup() %>% 
  mutate(size = as.integer(stringr::str_extract(size_class, "\\d+")))
#> drop_na: no rows removed
#> drop_na: removed 32 rows (2%), 1,920 rows remaining
#> group_by: 5 grouping variables (id, year, quarter, area, size_class)
#> summarise: now 120 rows and 6 columns, 4 group variables remaining (id, year, quarter, area)
#> ungroup: no grouping variables
#> mutate: new variable 'size' (integer) with 8 unique values and 0% NA

ggplot(levin_cod, aes(size, size_class)) + geom_point()

# Plot Levin
p1 <- levin_cod %>% filter(size > 5 & size < 60) %>% 
  ggplot(., aes(size, levin)) +
  geom_point() +
  stat_smooth() + 
  facet_wrap(~ area, scales = "free", ncol = 1) + 
  ggtitle("Cod")
#> filter: removed 46 rows (18%), 206 rows remaining
  
p2 <- levin_fle %>% filter(size > 5 & size < 40) %>% 
  ggplot(., aes(size, levin)) +
  geom_point() +
  stat_smooth() + 
  facet_wrap(~ area, scales = "free", ncol = 1) +
  ggtitle("Fle")
#> filter: removed 3 rows (2%), 117 rows remaining

p1 + p2
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'

# How can I know this is not simply due to samples size? (more samples, more diversity)
# Plot n prey in individual stomachs vs sample size in the year, quarter, ices_rect and size_class
# If there's a strong positive relationship, larger mean number of preys 
saturation_fle_prey <- fle_prey %>%
  mutate(size_class = cut(pred_cm, breaks = c(seq(0, 100, 5)))) %>% 
  pivot_longer(15:30) %>% 
  filter(value > 0) %>% 
  group_by(unique_pred_id, year, quarter, area, size_class) %>%
  summarise(n_prey = n()) %>% 
  ungroup() %>% 
  group_by(year, quarter, area, size_class) %>%
  summarise(mean_n_prey = mean(n_prey), 
            sample_size_n = n()) %>% 
  ungroup()
#> mutate: new variable 'size_class' (factor) with 8 unique values and 0% NA
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 2180x33, now 34880x19]
#> filter: removed 30,910 rows (89%), 3,970 rows remaining
#> group_by: 5 grouping variables (unique_pred_id, year, quarter, area, size_class)
#> summarise: now 1,785 rows and 6 columns, 4 group variables remaining (unique_pred_id, year, quarter, area)
#> ungroup: no grouping variables
#> group_by: 4 grouping variables (year, quarter, area, size_class)
#> summarise: now 120 rows and 6 columns, 3 group variables remaining (year, quarter, area)
#> ungroup: no grouping variables

ggplot(saturation_fle_prey, aes(sample_size_n, mean_n_prey, color = area)) +
  geom_jitter(alpha = 0.8, size = 2) +
  stat_smooth(method = "lm", se = FALSE) + 
  facet_wrap(~ size_class, scales = "free") + 
  scale_color_viridis(discrete = TRUE)
#> `geom_smooth()` using formula 'y ~ x'

saturation_cod_prey <- cod_prey %>%
  mutate(size_class = cut(pred_cm, breaks = c(seq(0, 100, 5)))) %>% 
  pivot_longer(15:30) %>% 
  filter(value > 0) %>% 
  group_by(unique_pred_id, year, quarter, area, size_class) %>%
  summarise(n_prey = n()) %>% 
  ungroup() %>% 
  group_by(year, quarter, area, size_class) %>%
  summarise(mean_n_prey = mean(n_prey), 
            sample_size_n = n()) %>% 
  ungroup()
#> mutate: new variable 'size_class' (factor) with 16 unique values and 0% NA
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x33, now 57488x19]
#> filter: removed 51,698 rows (90%), 5,790 rows remaining
#> group_by: 5 grouping variables (unique_pred_id, year, quarter, area, size_class)
#> summarise: now 3,229 rows and 6 columns, 4 group variables remaining (unique_pred_id, year, quarter, area)
#> ungroup: no grouping variables
#> group_by: 4 grouping variables (year, quarter, area, size_class)
#> summarise: now 252 rows and 6 columns, 3 group variables remaining (year, quarter, area)
#> ungroup: no grouping variables

saturation_cod_prey %>% 
  mutate(size = as.integer(stringr::str_extract(size_class, "\\d+"))) %>% 
  ggplot(., aes(size, mean_n_prey)) +
  geom_point() + 
  stat_smooth()
#> mutate: new variable 'size' (integer) with 16 unique values and 0% NA
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'

ggplot(saturation_cod_prey, aes(sample_size_n, mean_n_prey, color = area)) +
  geom_jitter(alpha = 0.8) +
  facet_wrap(~ size_class, scales = "free") + 
  stat_smooth(method = "lm", se = FALSE) + 
  scale_color_viridis(discrete = TRUE) +
  guides(color = FALSE)
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
#> `geom_smooth()` using formula 'y ~ x'